GROUP 1. Run title and other preliminaries
TEXT(Rayleigh-Taylor Instability      
TITLE
  DISPLAY
  A layer of liquid lies above a second layer of sligtly smaller
  density. The initial surface is not quite horizontal (effected
  by introducing a small gravity-related momentum source into the
  equation for U1), so motion begins.
 
  The representation is two-dimensional. The density is calculated
  from the enthalpy, by postulating a steep enthalpy-density
  dependence, modified by maximum and minimum limits.
                      _ _ _ _ _ _ _ _ _
                     |                 |
                     |      heavier    |
                  |  |..    liquid     |
                  |  |  :..        ..  |
                  |  |     :...  ..: :.|
                  |  |         :.:-----|----> interface
                  v  | lighter liquid  |
             gravity |_ _ _ _ _ _ _ _ _|
 
  The x-y plane is used, so that fields for all time steps can be
  dumped to a single phi file (called parphi) and viewed by PHOTON.
  The setting IDISPA=1 effects this. The PHOTON commands are
  supplied.
  ENDDIS
  photon use
  p
  parphi
 
 
  view z
  msg density contours
  do kk=1,20
  con rho1 z kk fi;0.1
  vec z kk
  enddo
  msg Density contours. Press return for another view
  pause;vec off
  view 1 1 1
  enduse
    GROUP 2. Transience; time-step specification
STEADY=F;GRDPWR(T,20,20.0,1.0)
    GROUP 3. X-direction grid specification
GRDPWR(X,20,1.0,1.0)
    GROUP 4. Y-direction grid specification
GRDPWR(Y,20,1.0,1.0)
    GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,V1,H1)
  ** Density is stored in order that it can be printed by way
     of the contour-plotting facility
STORE(RHO1)
    GROUP 8. Terms (in differential equations) & devices
TERMS(H1,N,Y,N,Y,Y,Y);DENPCO=T
    GROUP 9. Properties of the medium (or media)
  ** Set the temperature as TMP1A+TMP1B*H1
TMP1=LINH;TMP1A=0.0;TMP1B=1.0;cp1=1.0/tmp1b
  ** Density is made to depend steeply upon enthalpy; but the use of
     VARMAX and VARMIN prevents it from going outside narrow bounds.
     This practice preserves a sharp density difference, even though
     numerical diffusion spreads the enthalpy interface.
     Density is set to RHO1A+RHO1B*Temperature
RHO1=LINTEMP;RHO1A=50.0;RHO1B=-100.0
VARMAX(RHO1)=1.01;VARMIN(RHO1)=0.99
    GROUP 11. Initialization of variable or porosity fields
  ** lower layer
FIINIT(H1)=1.0;FIINIT(RHO1)=VARMIN(RHO1)
  ** upper layer
INIADD=F;PATCH(INIT,INIVAL,1,NX,NY/2+1,NY,1,1,1,1)
INIT(INIT,H1,0.0,0.0);INIT(INIT,RHO1,0.0,VARMAX(RHO1))
    GROUP 13. Boundary conditions and special sources
  ** Pressure relief
PATCH(REFP,CELL,1,1,1,1,1,1,1,LSTEP)
COVAL(REFP,P1,FIXP,0.0);COVAL(REFP,H1,ONLYMS,SAME)
COVAL(REFP,U1,ONLYMS,0.0);COVAL(REFP,V1,ONLYMS,0.0)
  ** Gravitational acceleration in v- and (to a lesser extent)
     u-directions.
PATCH(BUOYANCY,PHASEM,1,NX,1,NY,1,1,1,LSTEP)
COVAL(BUOYANCY,V1,FIXFLU,-9.81);COVAL(BUOYANCY,U1,FIXFLU,0.0981)
    GROUP 15. Termination of sweeps
LSWEEP=20;SELREF=T;RESFAC=0.01
    GROUP 21. Print-out of variables
  ** Print out pressure, velocity, temperature and density fields.
OUTPUT(P1,Y,Y,Y,Y,Y,Y);OUTPUT(U1,Y,Y,Y,Y,Y,Y)
OUTPUT(V1,Y,Y,Y,Y,Y,Y);OUTPUT(H1,Y,Y,Y,Y,Y,Y)
    GROUP 22. Spot-value print-out
IXMON=NX/4;IYMON=NY/2
SPEDAT(SET,GXMONI,TRANSIENT,L,F)
    GROUP 23. Field print-out and plot control
ITABL=1
PATCH(CONT,CONTUR,1,NX,1,NY,1,1,1,LSTEP)
PLOT(CONT,H1,0.0,10.0);PLOT(CONT,RHO1,0.0,20.0)
PLOT(CONT,P1,0.0,10.0);PLOT(CONT,U1,0.0,10.0)
PLOT(CONT,V1,0.0,10.0);ICHR=3
TSTSWP=-1;IDISPA=1