PHOTON USE
  p;;;;;
 
  norm
  gr ou x 1
  msg velocity vectors
  vec x 1 sh
  msg            -
  msg Press  to continue
  pause
  vec off;red
  msg pressure contours
  con p1 x 1 fi;0.005
  msg            -
  msg Press e to END
  enduse
 
    GROUP 1. Run title
TEXT(TURNAROUND DUCT WITH USER CHOICES  :B537
TITLE
  DISPLAY
  INVISCID OR TURBULENT compressible/incompressible
 
  2-dimensional (y-z), BFC, steady, elliptic simulation
 
               _____________________
             /                      --> outlet
            /     __________________
           |     |__________________              y ^
            \                       <-- inlet       |
              \ ____________________                +----> z
 
 
  This case concerns plane, two-dimensional, incompressible
  flow through a 180 degree turnaround duct. The flow can be
  treated as inviscid or turbulent, depending on user setting.
 
  enddis
 
REAL(WIN);WIN=1.0
    GROUP 6. Body-fitted coordinates or grid distortion
BFC=T;NONORT=T
  ** Create the body-fitted coordinate grid
ny=10;nz=20;xulast=50.0;yvlast=30.0;zwlast=100.0
GSET(D,nx,ny,nz,xulast,yvlast,zwlast)
  Shift plane k21 ysh1 from previous position in y direction
real(ysh1,ysh2,zsh,rotatex)
ysh1=110.0
GSET(C,K21,F,K21,+,0.0,ysh1,0.0)
  Shift plane k17 to zshift from k21 in z direction, including the
  intermediate points; power=1.0
zsh=-100.0
GSET(C,K17,F,K21,+,0.0,0.0,zsh,INC,1.0)
  Shift planes k5 with rotation in x direction of angle rotatex
  and in the y direction through the distance ysh2
rotatex=-3.14159;ysh2=95.0
GSET(C,K5,F,K17,RX,rotatex,ysh2,0.0,INC,1.0)
  Shift plane k1 -zsh in the z direction to complete grid
GSET(C,K1,F,K5,+,0.0,0.0,-zsh,INC,1.0)
   ** Set wup=t to account better for the high curvature of
      the w resolute...
WUP=T
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1);SOLUTN(P1,Y,Y,Y,N,N,N)
    GROUP 9. Properties of the medium (or media)
ENUL=0.0;ENUT=0.0
goto test
label loopback
    GROUP 11. Initialization of variable or porosity fields
FIINIT(P1)=1.E-10;FIINIT(W1)=WIN
 
    GROUP 13. Boundary conditions and special sources
patch(INLET,LOW,1,nx,1,NY,1,1,1,1)
coval(INLET,P1,fixflu,WIN);coval(inlet,W1,onlyms,WIN)
 
PATCH(OUTLET,HIGH,1,nx,1,NY,NZ,NZ,1,1);COVAL(OUTLET,P1,1.E4,0.0)
COVAL(OUTLET,V1,ONLYMS,0.0);COVAL(OUTLET,W1,ONLYMS,0.0)
 
    GROUP 15. Termination of sweeps
LSWEEP=200;SELREF=T; RESFAC=1.E-4
 
    GROUP 16. Termination of iterations
LITER(P1)=20
    GROUP 17. Under-relaxation devices
RELAX(P1,LINRLX,0.5);RELAX(V1,FALSDT,10.0);RELAX(W1,FALSDT,10.0)
mesg(Ask EXPERT to adjust relaxation parameters? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y) then
 #$030
endif
     GROUP 22. Spot-value print-out
IYMON=5;IZMON=10;NPRMON=LSWEEP
    GROUP 23. Field print-out and plot control
NPRINT=LSWEEP;ITABL=2;NPLT=5
PATCH(YX,CONTUR,1,1,1,NY,1,NZ,1,1)
PLOT(YX,P1,0.0,20.0)
PATCH(INNER,PROFIL,1,1,1,1,1,NZ,1,1)
PLOT(INNER,P1,0.0,0.0);PLOT(INNER,W1,0.0,0.0)
PATCH(OUTER,PROFIL,1,1,NY,NY,1,NZ,1,1)
PLOT(OUTER,P1,0.0,0.0);PLOT(OUTER,W1,0.0,0.0)
tstswp=-1;
goto stop
label test
char(vis)
 
+ do ii=1,5
+   mesg(
+ enddo
mesg( Initial data that can be changed:
mesga( Incompressible flow
mesg( Inviscid flow
mesg( Laminar flow
mesg( Inlet velocity is set to 1. m/s
mesga( Do you want to change settings (y/n)? (Default n)
readvdu(ans,char,n)
if(:ans:.eq.y) then
+ do ii=1,5
+   mesg(
+ enddo
+ mesg(Compressible flow ? (y/n)
+ readvdu(ans,char,n)
+ if(:ans:.eq.y) then
+   STORE(RHO1)
+   TERMS(W1,Y,Y,N,Y,Y,Y);TERMS(V1,Y,Y,N,Y,Y,Y)
  ** Set density to RHO1A*(P1+PRESS0)**RHO1B+RHO1C
+   RHO1=COMPRESS; RHO1A=1.0; RHO1B=0.714;PRESS0=1.0; RHO1C=0.0
  ** Set d(ln(density))/dp=1.0/(PRESS0/RHO1B+P1/RHO1B)
+   DRH1DP=COMPRESS
+   FIINIT(RHO1)=1.0
+ endif
 
+ do ii=1,5
+   mesg(
+ enddo
+ mesga(Introduce viscosity ? (y/n)
+ readvdu(vis,char,n)
+ if(:vis:.eq.y) then
+   enul=1.0
+   mesg(viscosity introduced and wall functions activated)
+   patch(innr,swall,1,nx,1,1,1,nz,1,1)
+   coval(innr,w1,1.0,0.0)
+   patch(outr,nwall,nx,1,ny,ny,1,nz,1,1)
+   coval(outr,w1,1.0,0.0)
+   char(turb)
+   mesg(Introduce turbulence? (y/n)
+   readvdu(turb,char,n)
+   if(:turb:.eq.y) then
+     turmod(kemodl)
+     fiinit(ke)=win**2*0.01
+     fiinit(ep)=fiinit(ke)**1.5/yvlast
+     patch(inlet,low,1,nx,1,ny,1,1,1,1)
+     coval(inlet,ke,onlyms,fiinit(ke))
+     coval(inlet,ep,onlyms,fiinit(ep))
+     patch(innr,swall,1,nx,1,1,1,nz,1,1)
+     coval(innr,w1,loglaw,0.0)
+     coval(innr,ke,loglaw,loglaw)
+     coval(innr,ep,loglaw,loglaw)
+     patch(outr,nwall,1,nx,ny,ny,1,nz,1,1)
+     coval(outr,w1,loglaw,0.0)
+     coval(outr,ke,loglaw,loglaw)
+     coval(outr,ep,loglaw,loglaw)
+   endif
+ endif
+ if(enul.gt.0.0) then
+   do ii=1,5
+     mesg(
+   enddo
+   mesga(Reynolds number = :win*yvlast/enul: OK?
+   mesg(If not, insert new value.
+   readvdu(WIN,real,1.0)
+   win=win*enul/yvlast
+ endif
endif
 
NZPRIN=NZ/5;NYPRIN=NY/5
WALPRN=T
 
goto loopback
label stop