GROUP 1. Run title and other preliminaries
TEXT(3D Steady Heat Conduction In Cube 
#cls
TITLE
  DISPLAY
                  /------/------/------/------/------/
                 /temperature fixed at 1.0 ----> X  /|
                /------/------/------/------/------/ /
               /      /      /      /      /      /|/|
              /------/------/------/------/------/ / /
             /      /      /      /      /      /|/|/|
            /------/------/------/------/------/ / / /
           /      /      /      /      /      /|/|/|/|
  north   |------|------|------|------|------| / / / /
    ^     |      |      |      |      |      |/|/|/|/|
    |     |------|------|------|------|------| / / / /  high
    |     |      |      |      |      |      |/|/|/|/    ^
    y     |------|------|------|------|------| / / /    /
    |     |      |      |      |      |      |/|/|/    z
    |     |------|------|------|------|------| / /    /
    |     |      |      |      |      |      |/|/    /
    |     |------|------|------|------|------| /    /
    |     |  X <-- temperature fixed at zero |/    /
  south   |------|------|------|------|------|   low
          /west----------------x----->east
#pause
    The following PHOTON "use" file enables the temperature
    distribution to be visualized. 
    
    PHOTON USE
    ext;;;;
    
    view 5 2 -1
    msg a cube with fixed temperatures at opposite corners
    msg press enter
    gr x m; gr  y m; gr  z 1
    pause;gr off;gr ou x m;gr ou y m;gr ou z 1;red
    msg the temperature for z=1
    msg press enter
    con temp z 1 fi;0.01
    con temp z 1 ;int 40                       
    gr ou x m;gr ou y m;gr ou z 1;
    pause;red

    msg the temperature for y=ny
    msg press enter
    con temp y m fi;0.01
    con temp y m ;int 40
    gr ou x m;gr ou y m;gr ou z 1;
    pause;red
    msg the temperature for x=nx
    msg press enter
    con temp x m fi;0.01
    con temp x m ;int 40
    gr ou x m;gr ou y m;gr ou z 1;
    pause;red
    msg Press e to END, or ? for help 
    enduse
  ENDDIS
#pause 
integer(lit)
REAL(XLENGTH,YLENGTH,ZLENGTH)
XLENGTH=1.0;YLENGTH=1.0;ZLENGTH=1.0
 
  domain size and grid
NX=10; NY=10; NZ=10; 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
IF(.NOT.:ANS:.EQ.Y) THEN
 SOLUTN(H1,Y,Y,Y,N,N,N);NAME(H1)=TEMP
ENDIF
   ** Store BLOK, the marker variable needed for the activation
      of the block-correction feature in the linear-equation
      solver.
STORE(BLOK)
   ** Indicate that it is to variable 14 that the block-
      correction feature will be applied
IVARBK=14
 
    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
 
   ** Define eight block-correction blocks by ascribing BLOK
     values to chosen locations by way of:-
     FIINIT, for block 1, and PATCH and INIT for the other
     blocks
     Note that INIADD is = F, so that the values set by INIT
     replace the value set by FIINIT. The only region which
     retains that value is that in the east, north, high
     corner.
 
INIADD=F;FIINIT(BLOK)=1.0
  West, south, low corner
PATCH(BLOK2,INIVAL,1,NX/2,1,NY/2,1,NZ/2,1,LSTEP)
INIT(BLOK2,BLOK,0.0,2.0)
  East, north, low corner
PATCH(BLOK3,INIVAL,NX/2+1,NX,1,NY/2,1,NZ/2,1,LSTEP)
INIT(BLOK3,BLOK,0.0,3.0)
  West, south, low corner
PATCH(BLOK4,INIVAL,1,NX/2,NY/2+1,NY,1,NZ/2,1,LSTEP)
INIT(BLOK4,BLOK,0.0,4.0)
  East, north, low corner
PATCH(BLOK5,INIVAL,NX/2+1,NX,NY/2+1,NY,1,NZ/2,1,LSTEP)
INIT(BLOK5,BLOK,0.0,5.0)
  West, south, high corner
PATCH(BLOK6,INIVAL,1,NX/2,1,NY/2,NZ/2+1,NZ,1,LSTEP)
INIT(BLOK6,BLOK,0.0,6.0)
  East, south, high corner
PATCH(BLOK7,INIVAL,NX/2+1,NX,1,NY/2,NZ/2+1,NZ,1,LSTEP)
INIT(BLOK7,BLOK,0.0,7.0)
  West, north, high corner
PATCH(BLOK8,INIVAL,1,NX/2,NY/2+1,NY,NZ/2+1,NZ,1,LSTEP)
INIT(BLOK8,BLOK,0.0,8.0)
 
    GROUP 13. Boundary conditions and special sources
   **Corner at IX=IY=IZ=1
PATCH(COLD,CELL,1,1,1,1,1,1,1,1)
   **Fix temperature to zero
COVAL(COLD,TEMP,1.E2,0.0)
   **Corner at IX=NX, IY=NY, IZ=NZ
PATCH(HOT,CELL,NX,NX,NY,NY,NZ,NZ,1,1)
   **Fix temperature to 1.0
COVAL(HOT,TEMP,1.E2,1.0)
 
    GROUP 15. Termination of sweeps
RESREF(TEMP)=1.E-6;LSWEEP=2
    GROUP 16. Termination of iterations
   **Terminate iterations when 500 iterations have been performed.
     (The minus sign is a signal to activate progress-of-solution
     print-out on the screen at a sweep frequency dictated by the
     value of TSTSWP.)
LITER(TEMP)=-500
  ** Set the over-relaxation factor for the linear-equation
     solver.
OVRRLX=1.7
  ** Set the frequency of application of the block=correction
     feature to once per iteration
ISOLBK=1
  ** 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

  ************************************************************
  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
PATCH(YLINE,PROFIL,NX/2,NX/2,1,NY,NZ/2,NZ/2,1,1)
PLOT(YLINE,TEMP,0.0,0.0);PLOT(YLINE,BLOK,0.0,0.0)
 
   **Plot a contour diagram for the plane IX=nx/2
PATCH(XPLANE,CONTUR,NX/2,NX/2,1,NY,1,NZ,1,1)
   **Let the diagram have 20 temperature intervals
PLOT(XPLANE,TEMP,0.0,20.0);PLOT(XPLANE,BLOK,0.0,20.0);ICHR=2
 
    GROUP 24. Dumps for restarts
 
       Interesting variants include the following:-
    1. By putting FIINIT(TEMP)= ... in GROUP 11, see how the
       initial guess influences solution time, final result etc.
    2. By changing NX, NY and NZ, see how increasing or decreasing
       the fineness of spatial subdivision influences solution
       time, results, etc.
    3. By changing XLENGTH, YLENGTH and ZLENGTH, see how the
       solution changes when the block is made long and thin,
       short and fat, etc.
    4. By changing PRNDTL(TEMP), see how the value of the
       conductivity changes the solution.
    5. By changing the third, fourth, and further  arguments of
       PATCH(HOT,... and PATCH(COLD,..., explore how changing the
       boundary-condition locations influence the solution.
       Also, introduce more PATCHes at which boundary conditions
       are provided; and display the results in different ways by
       introducing more PATCHes for profile and contour plots.
    6. By changing the COVAL coefficient from 1.E2 to FIXFLU,
       explore the effect of having fixed-flux rather than fixed-
       value boundary conditions; but remember that if all the
       boundary conditions are of fixed-flux type, the solution is
       not uniquely defined.
    7. By changing the values of OVRRLX, ISOLBK, ISOLX, ISOLY,
       ISOLZ, see how these solution-control parameters
       influence the solution and the computer time and number
       of iterations needed to attain it.
    8. By changing the number and location of the blocks used in
       the block-correction feature, establish the influences of
       factors on the accuracy and computer time.