TALK=F;RUN( 1, 1)
TEXT(USP. 3D Steady Heat Conduction In Cube
  created by dbs on 11.12.05 and modified by gin
  Its purpose is to compare the 'Stone-like' and usp solvers for a
  problem with a simple one-dimensional solution,
  A cube of conducting material contains a cube of non-conducting material.
  NX, NY and NZ can be varied.
  The case has been created by changing the boundary conditions of library 
  case 100, by inserting the non-conducting block, and by removing 
  distracting items such as BLOK.
  
   
REAL(XLENGTH,YLENGTH,ZLENGTH)
XLENGTH=1.0;YLENGTH=1.0;ZLENGTH=1.0
 
USP=T
UAUTO  = F
USPDBG=F
UTCPLT=F
USPVTK = T
USPIMB=F
MXLEV=0
MYLEV=0
MZLEV=0
DOMAT=-1

CELLST=10
FACEST=10
MINPRP=-1
MAXPRP=100

PARSOL = F

  domain size and grid
NX=20; NY=20; NZ=20
mesg(Default values of number of cells in X direction is :NX:
mesg(Do you want to change it? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
 mesg(Please, enter values of number of cells in X direction.
 readvdu(NX,int,NX)
endif

mesg(Default values of number of cells in Y direction is :NY:
mesg(Do you want to change it? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
 mesg(Please, enter values of number of cells in Y direction.
 readvdu(NY,int,NY)
endif

mesg(Default values of number of cells in Z direction is :NZ:
mesg(Do you want to change it? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
 mesg(Please, enter values of number of cells in Z direction.
 readvdu(NZ,int,NZ)
endif

BOOLEAN(lVert)
lVert=F
mesg(Do you want to view results in the centres of cells? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
lVert=T
endif

if(lVert)then
SPEDAT(SET,USPIO,VERTCENT,L,F)
endif


 xulast=xlength; yvlast=ylength; zwlast=zlength
#unigrid
label top
 
    GROUP 7. Variables stored, solved & named
   **Choose first-phase enthalpy (H1) as dependent variable
     and activate the whole-field elliptic solver
NAME(150)=TEMP
NAME(149)=PRPS
SOLVE(TEMP)
STORE(PRPS) 
FIINIT(PRPS)= -1
    GROUP 8. Terms (in differential equations) & devices
   **For pure conduction, cut out built-in source and convection
     terms
TERMS(TEMP,N,N,Y,N,Y,Y)
 
    GROUP 9. Properties of the medium (or media)
   **Thermal conductivity will be ENUL*RHO1/PRNDTL(TEMP), so :
ENUL=1.0;RHO1=1.0;PRNDTL(TEMP)=1.0
 
    GROUP 11. Initialization of variable or porosity fields
PATCH(CENTRAL, INIVAL, 2,NX-1,2,NY-1,1,NZ,1,1)
COVAL(CENTRAL, PRPS, 0.0, 199) 
    GROUP 13. Boundary conditions and special sources
   ** SLAB AT IZ=1
PATCH(COLD,CELL,1,NX,1,NY,1,1,1,1)
   **Fix temperature to zero
COVAL(COLD,TEMP,1.E5,0.0)
   ** SLAB AT IZ=NZ
PATCH(HOT,CELL,1,NX,1,NY,NZ,NZ,1,1)
   **Fix temperature to 1.0
COVAL(HOT,TEMP,1.E5,1.0)
 
    GROUP 15. Termination of sweeps
RESREF(TEMP)=1.E-6
ENDIT(TEMP)=1.e-6
    GROUP 16. Termination of iterations
LITER(TEMP)=1000
  ** Set the frequencies of application of the one-dimensional
     correction features in the linear-equation solver to once
     per iteration for each direction.
      ISOLX=1;ISOLY=1;ISOLZ=1
 RESFAC=1.E-7
  ************************************************************
  Group 17. Relaxation
RELAX(TEMP,LINRLX,1.0) 
 
    GROUP 21. Print-out of variables
   **Print fields of temperature
OUTPUT(TEMP,Y,N,N,N,N,N)
 
    GROUP 22. Spot-value print-out
IXMON=NX/2+1;IYMON=NY/2+1;IZMON=NZ/2+1;UWATCH=T
 
    GROUP 23. Field print-out and plot control
NXPRIN=NX/5;NYPRIN=NY/5;NZPRIN=NZ/5
   **Plot a profile along the line IX=nx/2,IZ=nz/2
 
   **Plot a contour diagram for the plane IX=nx/2
   **Let the diagram have 20 temperature intervals
ICHR=2
LSWEEP=100
NPRINT=1
YZPR=t
NXPRIN=NX/2
NZPRIN=1
NYPRIN=NY/2 
    GROUP 24. Dumps for restarts