TALK=F;RUN( 1, 1)

TEXT(2D HEAT CONDUCTION WITH SPACE DEPENDENT SOURCE
title
  DISPLAY 
    This case solves a two-dimensional steady diffusion problem
    with space-dependent source and boundary conditions.
    It can be tested against exact solution fields.

    Equation:
    d(dT/dX)/dX + d(dT/dX)/dX +  S(X,Y)  =  0
    S(X,Y)  =  2*sinX*sinY
    Boundary conditions:
    Y = 0 ;  T = 0
    Y = 1 ;  T = sin(1)*sin(X)
    X = 0 ;  T = 0
    X = 1 ;  T = sin(1)*sin(Y)
    Exact solution:
    T(X,Y) = sin(X)*sin(Y)    
    
  ENDDIS
  PHOTON USE
  p
 
 

 
  gr ou z 1
  msg           Temperature contours
  msg           (Numerical solution)
  con tem1 z 1 fi;0.001
  msg Press Enter to continue
  pause
  con off;red
  con exac z 1 fi;0.001
  msg           Temperature contours
  msg             (Exact solution)
  msg Press e to END
  ENDUSE
    GROUP 1. Run title and other preliminaries
REAL(XLENGTH,YLENGTH)
XLENGTH=1.0;YLENGTH=1.0
NX=25;NY=25
    GROUP 3. X-direction grid specification
GRDPWR(X,NX,XLENGTH,1.0)
    GROUP 4. Y-direction grid specification
GRDPWR(Y,NY,YLENGTH,1.0)
    GROUP 7. Variables stored, solved & named
SOLVE(TEM1); STORE(EXAC)
    GROUP 8. Terms (in differential equations) & devices
   For pure conduction, cut out built-in source and convection terms
TERMS(TEM1,N,N,Y,N,Y,Y)
    GROUP 9. Properties of the medium (or media)
PRNDTL(TEM1)=1;ENUL=1.
    GROUP 11. Initial values
REAL(XCOR,YCOR,DXC,DYC)
DXC=XULAST/NX
DYC=YVLAST/NY

  Following commands set the initial field as.. TEM1 = X + Y
DO II=1,NX
 XCOR=0.5*DXC+(II-1)*DXC
 DO JJ=1,NY
+ YCOR=0.5*DYC+(JJ-1)*DYC
+ PATCH(I:II:I:JJ:,INIVAL,:II:,:II:,:JJ:,:JJ:,1,NZ,1,1)
+ INIT(I:II:I:JJ:,EXAC,0.0,SIN(XCOR)*SIN(YCOR))
 ENDDO
ENDDO
    GROUP 13. Boundary conditions and special sources
   TEM1=0 at IX=0
PATCH(XEQ0,WWALL,1,1,1,NY,1,1,1,1)
   Fix temperature to zero
COVAL(XEQ0,TEM1,1./PRNDTL(TEM1),0.0)
   H1=0 at IY=0
PATCH(YEQ0,SWALL,1,NX,1,1,1,1,1,1)
   Fix temperature to zero
COVAL(YEQ0,TEM1,1./PRNDTL(TEM1),0.0)

   H1=sin1*sinX at Y=1
DO II=1,NX
 XCOR=0.5*DXC+(II-1)*DXC
 PATCH(YEQ1:II:,NWALL,:II:,:II:,NY,NY,1,1,1,1)
 COVAL(YEQ1:II:,TEM1,1./PRNDTL(TEM1),SIN(1.)*SIN(XCOR))
ENDDO
   H1=sin1*sinY at X=1
DO II=1,NY
 YCOR=0.5*DYC+(II-1)*DYC
 PATCH(XEQ1:II:,EWALL,NX,NX,:II:,:II:,1,1,1,1)
 COVAL(XEQ1:II:,TEM1,1./PRNDTL(TEM1),SIN(1.)*SIN(YCOR))
ENDDO

   Space dependent source
DO II=1,NX
 XCOR=0.5*DXC+(II-1)*DXC
 DO JJ=1,NY
+ YCOR=0.5*DYC+(JJ-1)*DYC
+ PATCH(S:II:S:JJ:,VOLUME,:II:,:II:,:JJ:,:JJ:,1,1,1,1)
+ COVAL(S:II:S:JJ:,TEM1,FIXFLU,2.*SIN(XCOR)*SIN(YCOR))
 ENDDO
ENDDO

    GROUP 15. Termination of sweeps
LSWEEP=10
    GROUP 21. Print-out of variables
OUTPUT(TEM1,Y,Y,Y,Y,Y,Y)
    GROUP 22. Spot-value print-out
IXMON=NX/2+1;IYMON=NY/2+1;IZMON=NZ/2+1
    GROUP 23. Field print-out and plot control
NXPRIN=NX/5;NYPRIN=NY/5;NZPRIN=NZ/5;nplt=1

  SAVE5BEGIN 
 (STORED of EXT1 is SIN(XG)*SIN(YG))
 (STORED of DTEM is EXT1-TEM1)  
  SAVE5END



      Usp related variables
USP    = T
UAUTO  = F
USPDBG = F
UTCPLT = F
USPVTK = T
USPIMB = F
MXLEV  = 4 
MYLEV  = 4 
MZLEV  = 4
DOMAT  = -1
MINPRP = -1
MAXPRP = 250
CELLST = 10
FACEST = 10
PARSOL = F
ISG62 = 0
ISG60 = 1
mesg(Do you want to view results in the centres of cells? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y)then
SPEDAT(SET,USPIO,VERTCENT,L,F)
endif

STOP