TALK=F;RUN( 1, 1)
TEXT(Flow in porous media. USP Test 6.                 
title
  DISPLAY
  
  This case solves a three-dimensional steady hydrodynamics
  problem with flow in propous medium.

  The analytical solution is: U1 = Uin, P1 = P-A*X
  where A depends on viscosity friction factor.
  ENDDIS  
real(fric,visc,Velin)
fric=10.
mesg(Default value of friction factor is :fric:
mesg(Do you want to change it? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
 mesg(Enter new value of friction factor
 readvdu(fric,real,fric)
 mesg(New value of friction factor is :fric:
endif

visc=0.
mesg(Default value of viscosity is :visc:
mesg(Do you want to change it? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
 mesg(Enter new value of viscosity 
 readvdu(visc,real,visc)
 mesg(New value of viscosity is :visc:
endif

Velin=1.
mesg(Default value of inlet velocity is :Velin: m/s
mesg(Do you want to change it? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
 mesg(Enter new value of velocity 
 readvdu(Velin,real,Velin)
 mesg(New value of velocity is :Velin: m/s
endif

STEADY=T

BOOLEAN(dim1)
dim1=F
mesg(Do you want to use 1D case (y) or 3D(n)? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
 dim1=T
endif

RSET(D,DOM,1.,1.,1.)
RHO1    = 1.
ENUL    = visc

IF(dim1)THEN
 RSET(M,10,1,1)
 STORE(V1,W1)
 SOLVE(P1,U1)
 SOLUTN(P1  ,Y,Y,Y,N,N,Y)
 SOLUTN(U1  ,Y,Y,Y,N,N,Y)
 LSWEEP = 5

 PATCH(IN,WEST,1,1,1,NY,1,NZ,1,1)
 COVAL(IN,P1, FIXFLU, Velin*RHO1)
 COVAL(IN,U1, 0., Velin)
 COVAL(IN,V1, 0., 0.)
 COVAL(IN,W1, 0., 0.)

 PATCH(OUT,EAST,NX,NX,1,NY,1,NZ,1,1)
 COVAL(OUT,P1, 1.E+3, 0.)

ELSE


RSET(M,10,10,10)

  Group 7. Variables: STOREd,SOLVEd,NAMEd
ONEPHS = T
SOLVE(P1,U1,V1,W1)
SOLUTN(P1  ,Y,Y,Y,N,N,Y)
SOLUTN(U1  ,Y,Y,Y,N,N,Y)
SOLUTN(V1  ,Y,Y,Y,N,N,Y)
SOLUTN(W1  ,Y,Y,Y,N,N,Y)

 
integer(idir)
idir = 0
mesg(Default direction of flow is from WEST to EAST
mesg(Do you want to change it? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
 mesg(Enter new direction of flow 
 mesg(0 = WEST-EAST, 1 = EAST-WEST
 mesg(2 = SOUTH-NORTH, 3 = NORTH-SOUTH
 mesg(4 = LOW-HIGH, 5 = HIGH-LOW
 readvdu(idir,int,0)
 if(idir.gt.5)then
 idir = 0
 endif
 if(idir.lt.0)then
 idir = 0
 endif
endif
         Flow from WEST to EAST
if(idir.eq.0)then
 mesg(Flow from WEST to EAST
 PATCH(IN,WEST,1,1,1,NY,1,NZ,1,1)
 COVAL(IN,P1, FIXFLU, Velin*RHO1)
 COVAL(IN,U1, 0., Velin)
 COVAL(IN,V1, 0., 0.)
 COVAL(IN,W1, 0., 0.)

 PATCH(OUT,EAST,NX,NX,1,NY,1,NZ,1,1)
 COVAL(OUT,P1, 1.E+3, 0.)
endif
         Flow from EAST to WEST
if(idir.eq.1)then
 mesg(Flow from EAST to WEST
 PATCH(IN,EAST,NX,NX,1,NY,1,NZ,1,1)
 COVAL(IN,P1, FIXFLU, -Velin*RHO1)
 COVAL(IN,U1, 0., -Velin)
 COVAL(IN,V1, 0., 0.)
 COVAL(IN,W1, 0., 0.)

 PATCH(OUT,WEST,1,1,1,NY,1,NZ,1,1)
 COVAL(OUT,P1, 1.E+3, 0.)
endif
         Flow from SOUTH to NORTH
if(idir.eq.2)then
 mesg(Flow from SOUTH to NORTH
 PATCH(IN,SOUTH,1,NX,1,1,1,NZ,1,1)
 COVAL(IN,P1, FIXFLU, Velin*RHO1)
 COVAL(IN,U1, 0., 0.)
 COVAL(IN,V1, 0., Velin)
 COVAL(IN,W1, 0., 0.)

 PATCH(OUT,NORTH,1,NX,NY,NY,1,NZ,1,1)
 COVAL(OUT,P1, 1.E+3, 0.)
endif
         Flow from NORTH to SOUTH
if(idir.eq.3)then
 mesg(Flow from NORTH to SOUTH
 PATCH(IN,NORTH,1,NX,NY,NY,1,NZ,1,1)
 COVAL(IN,P1, FIXFLU, -Velin*RHO1)
 COVAL(IN,U1, 0., 0.)
 COVAL(IN,V1, 0., -Velin)
 COVAL(IN,W1, 0., 0.)

 PATCH(OUT,SOUTH,1,NX,1,1,1,NZ,1,1)
 COVAL(OUT,P1, 1.E+3, 0.)
endif
         Flow from LOW to HIGH
if(idir.eq.4)then
 mesg(Flow from LOW to HIGH
 PATCH(IN,LOW,1,NX,1,NY,1,1,1,1)
 COVAL(IN,P1, FIXFLU, Velin*RHO1)
 COVAL(IN,U1, 0., 0.)
 COVAL(IN,V1, 0., 0.)
 COVAL(IN,W1, 0., Velin)

 PATCH(OUT,HIGH,1,NX,1,NY,NZ,NZ,1,1)
 COVAL(OUT,P1, 1.E+3, 0.)
endif
         Flow from HIGH to LOW
if(idir.eq.5)then
 mesg(Flow from HIGH to LOW
 PATCH(IN,HIGH,1,NX,1,NY,NZ,NZ,1,1)
 COVAL(IN,P1, FIXFLU, -Velin*RHO1)
 COVAL(IN,U1, 0., 0.)
 COVAL(IN,V1, 0., 0.)
 COVAL(IN,W1, 0., -Velin)

 PATCH(OUT,LOW,1,NX,1,NY,1,1,1,1)
 COVAL(OUT,P1, 1.E+3, 0.)
endif

  Group 15. Terminate Sweeps
 LSWEEP = 25
mesg(Default value of LSWEEP is :LSWEEP:
mesg(Do you want to change it? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
 mesg(Enter new value of LSWEEP 
 readvdu(LSWEEP,int,LSWEEP)
 mesg(New value of LSWEEP is :LSWEEP:
endif

ENDIF

    Group 11.Initialise Var/Porosity Fields
INIADD  =    F
FIINIT(U1) =0.0; FIINIT(V1) =0.0
FIINIT(W1) =0.0; FIINIT(P1) =0.0

PATCH(FRICT   ,VOLUME,1,NX,1,NY,1,NZ,1,1)
COVAL(FRICT   ,U1  , fric, 0.000000E+00)
COVAL(FRICT   ,V1  , fric, 0.000000E+00)
COVAL(FRICT   ,W1  , fric, 0.000000E+00)
 
RESFAC  = 1.E-06
 
  Group 16. Terminate Iterations
 LITER (P1)=200 ;LITER (U1)=200 ;LITER (V1)=200
 LITER (W1)=200
 ENDIT(P1) = 1.e-9; ENDIT(U1) = 1.e-9; ENDIT(V1) = 1.e-9
 ENDIT(W1) = 1.e-9

  Group 17. Relaxation
 RELAX(P1  ,LINRLX, 1.0)
 RELAX(U1  ,FALSDT, 1.E+8)
 RELAX(V1  ,FALSDT, 1.E+8)
 RELAX(W1  ,FALSDT, 1.E+8)

 ECHO    =    T
  Group 22. Monitor Print-Out
 IXMON=NX-1 ;IYMON=1 ;IZMON=1
 NPRMON=100000
 NPRMNT=1
 TSTSWP=-1
 
  Group 23.Field Print-Out & Plot Control
 NXPRIN=1; NYPRIN =1
 NPRINT  =    100000
 ISWPRF  =         1 ;ISWPRL =    100000

      Usp related variables
USP    = T
USPDBG = F
UAUTO  = F
UTCPLT = T
USPVTK = T
USPIMB = F
MXLEV  = 4 
MYLEV  = 4 
MZLEV  = 4
DOMAT  = -1
MINPRP = -1
MAXPRP = 250
CELLST = 1
FACEST = 1
USPREL = 0.7
PARSOL = F
ISG62 = 0
ISG60 = 1
mesg(Do you want to use collocated arrangement (y) or staggered one (n)? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
 SPEDAT(SET,USP,METHOD,I,1)
 LSWEEP = 100
endif

IF(dim1)THEN
mesg(Do you want to use the debug mode (y)? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
 USPDBG = T
 SPEDAT(SET,USPDBG,PRINTCOEF,L,T) 
 SPEDAT(SET,USPDBG,PRINTASSP,L,T) 
endif
ENDIF

STOP