define(PASS,c c c( c c=====================================================================72 c ----- $1-pass ----- c c This is an m4 macro definition which will construct a pass c in the $1 direction. It is invoked by, for example, c P_A_S_S( x, y ) (without the "_") c in which case, all instances of "$1" will be replaced with "x", and c all "$2" by "y". c c Why has this been done? In order to update any given row, PPM requires c information from up to two rows above and below the row of interest. c This information must be from the same computational time, however, c as the computation progresses, the large storage arrad DDD becomes c a mixture of old timestep information and new timestep information c as updated zone averages (rhonu, pnu, unu, etc) are stored in DDD. c Therefore, a small working scratch space of five rows (a slab) will c be created. Consider an X-pass, then the slab would be dimensioned c (1-nbdy:nx+nbdy,5): c c ---------------------------------------- c islab=5 | > c ---------------------------------------- c islab=4 | > c ---------------------------------------- c islab=3 | > c ---------------------------------------- c islab=2 | > c ---------------------------------------- c islab=1 | > c ---------------------------------------- c c initially, islab=1 is taken from ddd(*,*,-1) c initially, islab=2 is taken from ddd(*,*,0) c initially, islab=3 is taken from ddd(*,*,1) c initially, islab=4 is taken from ddd(*,*,2) c initially, islab=5 is taken from ddd(*,*,3) c c Pointers to the correct slab rows (slabbb, slabb, slab, slabt, slabtt) c (for 2rows below, 1 row below, the row, one row above, 2 rows above) c which would have the values (1,2,3,4,5). After the call to ppm and the c updated values stored in DDD, the counters are cylcled to (2,3,4,5,1) c and the required data from DDD is written into the row pointed to c by the new value of slabtt. c c The user should also note that the velocities switch roles. That is, c during an X-pass, u1 is the x-velocity and u2 is the y-velocity while c during a Y-pass, u1 is the y-velocity and u2 is the x-velocity. c c This reduces both the scratch space requirement and the amount of c data copying required. c c the xvelocity is element #3 in DDD c the yvelocity is element #4 in DDD kvx = 3 kvy = 4 do i = 1-nbdy, n$1+nbdy+1 xxl(i) = $1l(i) enddo c c dddslab(*,*,1) = density c dddslab(*,*,2) = total energy c dddslab(*,*,3) = velocity c dddslab(*,*,4) = velocity c dddslab(*,*,5) = pressure c dddslab(*,*,6) = Eulerian sound speed c c initialize first 4 rows in slab c c ((do 100) do islab = 1, 4 i$2 = islab-2 do i$1 = 1-nbdy, n$1+nbdy dddslab(i$1, islab,1) = ddd( 1, ix, iy) dddslab(i$1, islab,2) = ddd( 2, ix, iy) dddslab(i$1, islab,3) = ddd(kv$1, ix, iy) dddslab(i$1, islab,4) = ddd(kv$2, ix, iy) enddo enddo c (do 100)) c c ((do 400) do i$2 = 1, n$2 c find pointers to proper rows in the slab to where c the sow two zones below are stored (slabbb) one c row below (slabb), etc islabbb = 1+mod(i$2-1,5) islabb = 1+mod(i$2, 5) islab = 1+mod(i$2+1,5) islabt = 1+mod(i$2+2,5) islabtt = 1+mod(i$2+3,5) c c At this point, we are ALWAYS able to fill position slabtt from ddd l$2 = i$2+2 do i$1 = 1-nbdy, n$1+nbdy l$1 = i$1 dddslab(i$1, islabtt,1) = ddd( 1, lx, ly) dddslab(i$1, islabtt,2) = ddd( 2, lx, ly) dddslab(i$1, islabtt,3) = ddd(kv$1, lx, ly) dddslab(i$1, islabtt,4) = ddd(kv$2, lx, ly) enddo ip = 1-nbdy d2inv = four / ($2l(i$2+1)+$2l(i$2+2)-$2l(i$2-1)-$2l(i$2)) $1courmx = zero do ijk = ip,n$1+nbdy 6 enddo call eos( dddslab(ip,islabbb,1), dddslab(ip,islabbb,2), ^ dddslab(ip,islabbb,3), dddslab(ip,islabbb,4), ^ dddslab(ip,islabbb,5), dddslab(ip,islabbb,6), ^ gamma, (n$1+2*nbdy) ) call eos( dddslab(ip,islab,1), dddslab(ip,islab,2), ^ dddslab(ip,islab,3), dddslab(ip,islab,4), ^ dddslab(ip,islab,5), dddslab(ip,islab,6), ^ gamma, (n$1+2*nbdy) ) call eos( dddslab(ip,islabtt,1), dddslab(ip,islabtt,2), ^ dddslab(ip,islabtt,3), dddslab(ip,islabtt,4), ^ dddslab(ip,islabtt,5), dddslab(ip,islabtt,6), ^ gamma, (n$1+2*nbdy) ) call ifelse(ifdouble,1,d_, ) do_ppmlr0_2dc ( xxl(ip), ^ dddslab(ip,islabbb,1), dddslab(ip,islab,1), dddslab(ip,islabtt,1), ^ dddslab(ip,islabbb,5), dddslab(ip,islab,5), dddslab(ip,islabtt,5), ^ dddslab(ip,islabbb,3), dddslab(ip,islab,3), dddslab(ip,islabtt,3), ^ dddslab(ip,islabbb,6), dddslab(ip,islab,6), dddslab(ip,islabtt,6), ^ dddslab(ip,islabbb,4), dddslab(ip,islabb,4), dddslab(ip,islab,4), ^ dddslab(ip,islabt,4), dddslab(ip,islabtt,4), ^ dddslab(ip,islab,2), ^ rhonu(ip), etotnu(ip), u1nu(ip), u2nu(ip), ^ dt, smlrho, smallp, smallu, ^ smalle, $1courmx, n$1, nbdy) c print *,' i$2, $1courmx=', i$2, $1courmx c courmx = max ( courmx, $1courmx ) c do i$1 = 1, n$1 ddd( 1,ix,iy) = rhonu(i$1) ddd( 2,ix,iy) = etotnu(i$1) ddd(kv$1,ix,iy) = u1nu(i$1) ddd(kv$2,ix,iy) = u2nu(i$1) enddo c enddo c (do 400)) c c c c ) c) C############################################################### c======================================================================= c= = c= To compile and run: = c= = c======================================================================= c c m4 Driverlagremap2dg0c.EXA > driver.f c c FORTRAN90: c f90 -trapuv -n32 -extend_source driver.f ppmlib98a.n32.a -lfortran cio.c c or c f90 -trapuv -64 -extend_source driver.f ppmlib98a.64.a -lfortran cio.c c c FORTRAN77: c f77 -trapuv -n32 -extend_source driver.f ppmlib98a.n32.a -lfortran cio.c c or c f77 -trapuv -64 -extend_source driver.f ppmlib98a.64.a -lfortran cio.c c c c======================================================================= define(IQX,256) define(IQY,256) define(IQ,256) define(PREFIX,L) define(zero,ifelse(ifdouble,0,0.,0.0d+00)) define(one,ifelse(ifdouble,0,1.,1.0d+00)) define(two,ifelse(ifdouble,0,2.,2.0d+00)) define(three,ifelse(ifdouble,0,3.,3.0d+00)) define(four,ifelse(ifdouble,0,4.,4.0d+00)) define(five,ifelse(ifdouble,0,5.,5.0d+00)) define(ten,ifelse(ifdouble,0,10.,1.0d+01)) define(half,ifelse(ifdouble,0,.5,5.0d-01)) define(qtr,ifelse(ifdouble,0,.25,2.5d-01)) define(eighth,ifelse(ifdouble,0,.125,1.25d-01)) define(fifth,ifelse(ifdouble,0,.2,2.0d-01)) c c The above macro definitions allow for double precision c======================================================================= c c THIS CODE PACKAGE IS SET UP ASSUMING THAT IT WILL FIRST BE c PREPROCESSED USING THE M4 UNIX MACRO PROCESSOR, BASED UPON c VALUES SET FOR THE ABOVE FLAGS AND CONSTANTS. c c======================================================================= c c c The constant IQX gives the number of zones of a physical row c (x-direction) (not counting the fake zones used to implement c boundary conditions.) c c The constant IQY gives the number of zones of a physical column c (y-direction) (not counting the fake zones used to implement c boundary conditions.) c c The constant IQ is the maximum of IQX, IQY, and For a given c problem, it is generally set to the maximum of the physical row c lengths in the x- and the y-directions. If IQX, IQY and IQ c are inconsistant, an error will result and the job will be terminated. c c PREFIX is a one character alpha-numeric. The names of all output image c dumps will be begin with this character. c c ************************************************************************** c * THIS CODE IS Driverlagremap2dg0c * c * THIS CODE WAS WRITTEN BY: WOODWARD RESEARCH GROUP. (c) 1999 * c * ALL RIGHTS RESERVED. * c ************************************************************************** c c c C WHAT This program is an example of a simple driver to c compute hydrodynamics in a lagrangian step + remap style c using routines from the PPMLIB library. Directional c splitting is used to approximate the 2-dimensional c Euler equations by applying a 1-dimensional operator c to treazt gradients in the x- and y- directions c separately. The inital problem consists of a two- c dimensional shock tube (discontinuities in both x- and c y-directions dividing the initial state into four c quadrants each with a different but constant initial c state. Reflecting boundary conditions are used. No c input is required. Just compile and run, linking with c the PPMLIB library. c c c ASSUMPTIONS 2 dimensions c c Cartesian coordinates c c Uniform grid c c No external body forces (e.g. no gravity) c c Gamma law equation of state: c p = (gamma - 1) * rho * ei c where c p = pressure c rho = density c ei = internal energy c gamma = ratio of specific heats, here gamma = 1.4 c c Initial state: c p(left) = 1.0; p(right) = 0.1; c rho(left) = 1.0; rho(right) = 0.125; c gamma = EOSGAM c initial velocity is everywhere zero. c initial discontinuity at center of grid. c c The real portion of the grid ranges from x=0.0 to c x=1.0 (XMAX). c c a fixed timestep is used. c The program runs for a fixed number of timesteps. C C BOUNDARIES The boundary conditions (here, reflecting) are handeled c by appending 9 "fake" zones to each end of the data c arrays. These zones are then filled with appropriate c values during a "boudary condition" step. c C C C DIMENSIONS c 1D: Two types of 1-d arrays exist: zone interface centered and c zone centered. An example of the first type would be an c array containing locations of the lefthand zone interfaces c XL(I). An example of the second would be zone widths, c DX(I) (=XL(I+1)-XL(I)), zone averaged density, pressure, c velocity, etc. Most arrays are zone centered. c c It will be assumed that the dimensions of all zone-averaged c type quantities are C [1-NBDY:N+NBDY], c where N is the number of real zones in the calculation, c NBDY is the number fake zones necessary for c boundary conditions. C The dimensions of all interface type quantities are, C [1-NBDY:N+NBDY+1]. C Note the "+1" please. c c 2D: The two-dimensional arrays containing the zone-averaged c density (rrho), pressure (pp), and velocities (uux, uuy) c have the dimensions c [1-NBDY:NX+NBDY, 1-NBDY:NY+NBDY]. c Note that the names of 2d arrays are "stuttered" (e.g. rrho) c From these are extracted working 1d arrays. c C output At the completion of each timestep, new zone-averaged C quantities [e.g. RHONU, the new zone averaged densities] c will have updated values on the range C [1:N], and will be stored in the c appropriate locations of the 2d arrays. C Any returned interface type thingies will have valid C values on the range C [1:N+1]. c c C OUTPUT QUANTITIES c output consists binary files scaled so that the values c are in the range 0-255. That is one byte per zone. c Each file is NXbyNY bytes giving a one pixel per c zone ratio when displayed with the XRAZ routines c available from LCSE. c c C program Driverlagremap2dg0c c c parameter (NXZONES = IQX) parameter (NYZONES = IQY) parameter (MAXSTEPS = 500) parameter (EOSGAM = 1.4) parameter (DT = 0.001250) parameter (XMAX = 1.0) parameter (YMAX = 1.0) parameter (N2DUMP = 20) parameter (T2DUMP = 99999.) c Note that because of the fixed timestep, N, XMAX, and DT are not c independent. c c NXZONES = number of real computational zones in the x-direction, set c by an m4 definition at beginning of code. c NYZONES = number of real computational zones in the y-direction, set c by an m4 definition at beginning of code. c NXZONES and NYZONES should be changed by changing the values of c IQX, IQY, and IQ at beginning. c MAXSTEPS = maximum number of timpsteps to perform. c EOSGAM = the equation of state gamma c DT = the fixed timestep c XMAX = the righthand edge x-value of the grid. c The x-grid runs from (0,XMAX) c YMAX = the righthand edge y-value of the grid. c The y-grid runs from (0,YMAX) c N2DUMP is the frequency of image dumps in iterations c T2DUMP is the frequency of image dumps in program time. c c Do not change this one: parameter (NBDY = 9) c c c======================================================================= c c The constant IQX gives the maximum length of a row in the c x-direction which will be encountered. c The constant IQY gives the maximum length of a column in the c y-direction which will be encountered. c c c======================================================================= c c ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension xl( 1-NBDY:NXZONES+NBDY+1) dimension yl( 1-NBDY:NYZONES+NBDY+1) dimension ddd( 1:4, 1-NBDY:NXZONES+NBDY, 1-NBDY:NYZONES+NBDY ) c c c c ddd = the 3-D array into which the 4 principal 2-D arrays of c basic physical variables are interleaved in order to c optimize memory bandwidth on microprocessor systems. c c xl = the locations of the lefthand zone interfaces c yl = the locations of the bottom zone interfaces c ddd = the 3-D array into which the 4 principal 2-D arrays of c basic physical variables are interleaved in order to c optimize memory bandwidth on microprocessor systems. c ddd(1,i,j) = rrho(i,j) c ddd(2,i,j) = eetot(i,j) c ddd(3,i,j) = uux(i,j) c ddd(4,i,j) = uuy(i,j) c rrho = the zone-averaged densities c eetot = the zone-averaged total energies c uux = the zone-averaged x-velocities c uuy = the zone-averaged y-velocities c c c smlrho = a trivial value of density, used to protect divides. c smalle = a trivial value of energy, used to protect divides. c smallp = a trivial value of pressure, used to protect divides. c smallu = a trivial value of velocity, used to protect divides. c smlrho = ifelse(ifdouble,0,1.0e-06, 1.0e-08) smalle = ifelse(ifdouble,0,1.0e-06, 1.0e-08) smallp = ifelse(ifdouble,0,1.0e-06, 1.0e-08) smallu = ifelse(ifdouble,0,1.0e-06, 1.0e-08) c c c c initialize the grid: xlimit = XMAX ylimit = YMAX c call setgrid( xl, yl, xlimit, ylimit, nxzones, nyzones, nbdy ) c c c initialize the hydrodynamic quantities: c call set2dshock( ddd, nxzones, nyzones, nbdy ) c c c nstop = MAXSTEPS tstop = float (MAXSTEPS) * DT ndump = N2DUMP tdump = T2DUMP c nstop is the maximum number of iterations to take for the entire c calculation c tstop is the maximum program time for the entire calculation, here c assuming a constant timestep DT c ndump is the frequency of image dumps in iteratiopns c tdump is the frequency of image dumps in program time. c c CONTROL performs the entire hydrodynamic calculation. c dtime = DT gamma = EOSGAM call cntrol ( xl, yl, ddd, gamma, dtime, smlrho, smalle, smallp, ^ smallu, ndump, tdump, nstop, tstop, nxzones, nyzones, nbdy) stop end subroutine cntrol ( xl, yl, ddd, gamma, dtime, smlrho, smalle, ^ smallp, smallu, ndump, tdump, nstop, tstop, nx, ny, nbdy) c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** ***** c ***** cc n n ttttt rrr ooo l ***** c ***** c c nn n t r r o o l ***** c ***** c n n n t rrr o o l ***** c ***** c n n n t rr o o l ***** c ***** c c n nn t r r o o l ***** c ***** cc n n t r r ooo llll ***** c ***** ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension xl( 1-NBDY:NX+NBDY+1) dimension yl( 1-NBDY:NY+NBDY+1) dimension ddd( 1:4, 1-NBDY:NX+NBDY, 1-NBDY:NY+NBDY ) c c c This routine runs the entire hydro calculation. c c c c ddd = the 3-D array into which the 4 principal 2-D arrays of c basic physical variables are interleaved in order to c optimize memory bandwidth on microprocessor systems. c nx = the number of real zones in the x-dimension for this c tile of the computational domain. c ny = the number of real zones in the y-dimension for this c tile of the computational domain. c nbdy = the number of fake zones needed to implement boundary c conditions so that in the hydrodynamical calculation c all zones can be treated as if they were interior zones. c c c dtime = the time step. c courmx = the maximum value of the Courant number for this tile c of the computational domain. c gamma = the value of the ratio of specific heats, as in the c gamma-law equation of state: p = (gamma-1)*rho*ei c c smlrho = a trivial value of density, used to protect divides. c smalle = a trivial value of energy, used to protect divides. c smallp = a trivial value of pressure, used to protect divides. c smallu = a trivial value of velocity, used to protect divides. c c xl = the locations of the lefthand zone interfaces c yl = the locations of the bottom zone interfaces c ddd = the 3-D array into which the 4 principal 2-D arrays of c basic physical variables are interleaved in order to c optimize memory bandwidth on microprocessor systems. c ddd(1,i,j) = rrho(i,j) c ddd(2,i,j) = eetot(i,j) c ddd(3,i,j) = uux(i,j) c ddd(4,i,j) = uuy(i,j) c rrho = the zone-averaged densities for this tile. c eetot = the zone-averaged total energies for this tile. c uux = the zone-averaged x-velocities for this tile. c uuy = the zone-averaged y-velocities for this tile. c c c Create an image of the initial state. The initial state will be c sequence number 0. numdmp = 0 call dump (ddd, numdmp, nx, ny, nbdy) c numdmp = the sequential number of the image 0->whatever, numdmp c will be automatically incremented after write the image. c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c The 1000 loop advances the calculation by two timesteps per c c iteration c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c initialize counters c ncycle = 0 time = 0.0 ntodmp = ndump ttodump = tdump c 1000 continue c c clear the old maximum courant number in preparation of the next pair of c timesteps: courmx = zero c c RUNHYD performs a symmetrized pair of passes x-y-y-x which c equals two complete timesteps. c call runhyd_lr ( xl, yl, ddd, dtime, courmx, gamma, 1 smlrho, smalle, smallp, smallu, 2 nx, ny, nbdy) c c c c update counters c ncycle = ncycle + 2 ntodmp = ntodmp - 2 time = time + dtime + dtime c write (6,1710) ncycle, time, courmx, dtime c c Check to see if time to write am image file. An image is c written every NDUMP steps or at itervals of DTDUMP in c problem time. c if (ntodmp) 1600,1600,1500 1500 if (time - tdump) 1700,1600,1600 1600 ntodmp = ndump ttodump = ttodump + tdump call dump (ddd, numdmp, nx, ny, nbdy) c numdmp = the sequential number of the image 0->whatever, numdmp c will be automatically incremented after write the image. call flush (6) c 1700 continue 1710 format (1h ,3hn =,i5,4x,3ht =, 1 1pe12.5,4x,13hcourant no. =,e12.5,4x,4hdt =,e12.5) c c Check to see if time to stop. Progam execution ends when either c the maximum number of timesteps, NSTOP, has been computed or c when the maximum problem time, TSTOP, has been reached. c if (ncycle .ge. nstop) stop if (time .ge. tstop) stop goto 1000 end c cSubroutineSubroutineSubroutineSubroutineSubroutineSubroutineSubroutineSubroutinc c #### # # ##### ##### #### # # ##### # # # ###### c # # # # # # # # # # # # # ## # # c #### # # ##### # # # # # # # # # # # ##### c # # # # # ##### # # # # # # # # # # c # # # # # # # # # # # # # # # ## # c #### #### ##### # # #### #### # # # # ###### c cSubroutineSubroutineSubroutineSubroutineSubroutineSubroutineSubroutineSubroutinc c subroutine setgrid( xl, yl, xmax, ymax, nx, ny, nbdy ) c c c c This routine sets the initial values of all necessary variables c in the vector common block. c This routine initializes the locations of the lefthand zone interfaces in c both the x- and y-directions, assuming a uniform grid. The grid consists c of NX zones (defined by NX+1 lefthand zone interfaces) in the x-direction c ranging from 0.0 to XMAX, with NBDY zones added to each end to allow for c boundary conditions. there are NY zones (defined by NY+1 lefthand zone c interfaces) in the y-direction ranging from 0.0 to YMAX, with NBDY zones c added to each end to allow for boundary conditions. c c c xl = the x-locations of the lefthand zone interfaces. c yl = the y-locations of the lefthand zone interfaces. c c xmax = The real (non boundary) portion of the grid is c XXL(1)=0.0 to XXL(NX+1)=XMAX. c ymax = The real (non boundary) portion of the grid is c YYL(1)=0.0 to YYL(NY+1)=YMAX. c c nx = number of real zones in in the x-direction c ny = number of real zones in in the y-direction c nbdy = number of boundary zones added to each end of the arrays to c implement boundary conditions. ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension xl(1-nbdy:nx+nbdy+1) dimension yl(1-nbdy:ny+nbdy+1) c deex = xmax / float(NX) DO I=1,NX+1 xl(i) = deex*float(i-1) ENDDO deey = ymax / float(NY) DO I=1,NY+1 yl(i) = deey*float(i-1) ENDDO return end subroutine set2dshock( ddd, nx, ny, nbdy ) c c c c This routine sets the initial values of all hydrodynamic c quantities for BKE's 2d shock tube problem. The initial c state consists of four different steady states: c c rhoq2 = 2.0 rhoq1 = 1.0 c pq2 = 2.0 pq1 = 1.0 c uxq2 = 0.0 uxq1 = 0.0 c uyq2 = 0.0 uyq1 = 0.0 c c rhoq3 = 0.5 rhoq4 = 0.1 c pq3 = 0.5 pq4 = 0.1 c uxq3 = 0.0 uxq4 = 0.0 c uyq3 = 0.0 uyq4 = 0.0 c c assuming a gamma law equation of state: p = (gamma-1)*rho*ei c c ..................................................... c . . . c . . . c . . . c . Q2 . Q1 . c . . . c . . . c ..................................................... < nydisc c . . . c . . . c . . . c . . . c . Q3 . Q4 . c . . . c . . . c . . . c . . . c ..................................................... c ^ c nxdisc c c ddd = the 3-D array into which the 4 principal 2-D arrays of c basic physical variables are interleaved in order to c optimize memory bandwidth on microprocessor systems. c ddd(1,i,j) = rrho(i,j) c ddd(2,i,j) = eetot(i,j) c ddd(3,i,j) = uux(i,j) c ddd(4,i,j) = uuy(i,j) c ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension ddd( 1:4, 1-NBDY:NX+NBDY, 1-NBDY:NY+NBDY ) c parameter( rhoq1 = 1.0) parameter( pq1 = 1.0) parameter( uxq1 = 0.0) parameter( uyq1 = 0.0) parameter( rhoq2 = 2.0) parameter( pq2 = 2.0) parameter( uxq2 = 0.0) parameter( uyq2 = 0.0) parameter( rhoq3 = 0.5) parameter( pq3 = 0.5) parameter( uxq3 = 0.0) parameter( uyq3 = 0.0) parameter( rhoq4 = 0.1) parameter( pq4 = 0.1) parameter( uxq4 = 0.0) parameter( uyq4 = 0.0) gammaq1 = 1.4 gammaq2 = 1.4 gammaq3 = 1.4 gammaq4 = 1.4 etq1 = 0.5*(uxq1**2+uyq1**2) + pq1/( (gammaq1-1.0)*rhoq1) etq2 = 0.5*(uxq2**2+uyq2**2) + pq2/( (gammaq2-1.0)*rhoq2) etq3 = 0.5*(uxq3**2+uyq3**2) + pq3/( (gammaq3-1.0)*rhoq3) etq4 = 0.5*(uxq4**2+uyq4**2) + pq4/( (gammaq4-1.0)*rhoq4) c c Location of discontinuity in the x-direction: nxdisc = nx/4 c Location of discontinuity in the y-direction: nydisc = (4*ny)/5 do 100 k = nydisc,ny+nbdy do 100 i = nxdisc,nx+nbdy ddd(1,i,k) = rhoq1 ddd(2,i,k) = etq1 ddd(3,i,k) = uxq1 ddd(4,i,k) = uyq1 100 continue c do 200 k = nydisc,ny+nbdy do 200 i = -nbdy+1,nxdisc-1 ddd(1,i,k) = rhoq2 ddd(2,i,k) = etq2 ddd(3,i,k) = uxq2 ddd(4,i,k) = uyq2 200 continue c do 300 k = 1-nbdy,nydisc-1 do 300 i = -nbdy+1,nxdisc-1 ddd(1,i,k) = rhoq3 ddd(2,i,k) = etq3 ddd(3,i,k) = uxq3 ddd(4,i,k) = uyq3 300 continue c do 400 k = 1-nbdy,nydisc-1 do 400 i = nxdisc,nx+nbdy ddd(1,i,k) = rhoq4 ddd(2,i,k) = etq4 ddd(3,i,k) = uxq4 ddd(4,i,k) = uyq4 400 continue c c c c return end subroutine dump (ddd, numdmp, nx, ny, nbdy ) c c ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) c c we are counting on the F77 compiler to automatically allocate/deallocate dimension realimg1(nx*ny) dimension realimg2(nx*ny) dimension realimg3(nx*ny) dimension realimg4(nx*ny) c c c ddd = the 3-D array into which the 4 principal 2-D arrays of c basic physical variables are interleaved in order to c optimize memory bandwidth on microprocessor systems. c ddd(1,i,j) = rrho(i,j) c ddd(2,i,j) = eetot(i,j) c ddd(3,i,j) = uux(i,j) c ddd(4,i,j) = uuy(i,j) c rrho = the zone-averaged densities. c eetot = the zone-averaged total energies c uux = the zone-averaged x-velocities. c uuy = the zone-averaged y-velocities. c c numdmp = the sequential number of the image 0->whatever, numdmp c will be automatically incremented after write the image. c nx = number of real computational zones in the x-direction, set c by an m4 definition at beginning of code. c ny = number of real computational zones in the y-direction, set c by an m4 definition at beginning of code. c c nbdy = the number fake zones necessary for boundary conditions. c dimension ddd( 1:4, 1-NBDY:NX+NBDY, 1-NBDY:NY+NBDY ) c c c c extract the non boundary zone data from DDD, placing them c contiguously into 1d arrays c ii = 0 do 1000 j = 1,ny do 1000 i = 1,nx ii = ii + 1 realimg1(ii) = ddd(1,i,j) realimg2(ii) = ddd(2,i,j) realimg3(ii) = ddd(3,i,j) realimg4(ii) = ddd(4,i,j) 1000 continue c c create a scale binary image file: c call wrtimg4 (realimg1, realimg2, realimg3, realimg4, 'PREFIX', 'r', ^ nx, ny, numdmp) numdmp = numdmp + 1 c c return end c c======================================================================= c c======================================================================= c subroutine wrtimg (t, acharc, acharc2, nx, ny, iimage, tmpmin, tmpmax) c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** ***** c ***** w w rrrr ttttt iii m m ggg ***** c ***** w w w r r t i m m g g ***** c ***** ww ww rrrr t i mm mm g ***** c ***** ww ww r r t i mm mm g gg ***** c ***** w w r r t i m m m g g ***** c ***** w w r r t iii m m m ggg ***** c ***** ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension t(nx,ny) parameter (nngx=IQX) parameter (nngy=IQY) parameter (ng=nngx*nngy) dimension igray(ng) character*8 namdmp character*1 acharc character*1 acharc2 c c t = the 2-D temparay array containing the data to write to disk c nx = the x dimension of the array t. c ny = the y dimension of the array t. c iimage = the sequential number of the image 0->whatever c igray = the 2-D array of grayscale values (ranging from 0 to 255) c acharc = the first character in the names of all images files c acharc2 = the second character in the names of all images files c c all binary image files are named 'abim####' where a = acharc, c b = acharc2, and #### is the sequential number 0000 to 9999 c namdmp = the 8-character name of the image file to be created. c c c Construct name of image file and store in NAMDMP: c write (namdmp,1004) acharc, acharc2, iimage 1004 format (a1,a1,'im',i4.4) c c c Scale array T so that it ranges from 0 to 255, storing the result in c the integer array IGRAY. Also returns the minimum (TMPMIN) and maximum c (TMPMAX) of array T. c call makgry (t, igray, tmpmin, tmpmax, nx, ny) c c Now write it to disk. c call wrtgry (igray, nx, ny, namdmp) c return end subroutine makgry (gray,igray,grymin,grymax,ngx,ngy) c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** ***** c ***** m m a k k ggg rrrr y y ***** c ***** m m a k k g g r r y y ***** c ***** mm mm a a kk g rrrr y y ***** c ***** mm mm a a k k g gg r r y ***** c ***** m m m aaaaa k k g g r r y ***** c ***** m m m a a k k ggg r r y ***** c ***** ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension gray(ngx,ngy), igray(ngx,ngy) c c This routine scales the data in the 2-D array GRAY and constructs c an array, IGRAY, of 2-byte grayscale levels. c c gray = the 2-D array containing the unscaled grayscale values. c igray = the 2-D array of grayscale values (ranging from 0 to 255) c ngx = the number of columns in the array gray. c ngy = the number of rows in the array gray. c c c Find the maximum and minimum: grymax = gray(1,1) grymin = gray(1,1) do j=1,ngy do i=1,ngx grymax = max (grymax,gray(i,j)) grymin = min (grymin,gray(i,j)) enddo enddo c c Now scale the array c denom = (grymax - grymin) if( denom .ne. 0.0 ) denom = 1.0/denom factor = 255.0 * denom c do 2000 j = 1,ngy do 1000 i = 1,ngx thyng = factor * (gray(i,j) - grymin) igray(i,j) = thyng 1000 continue 2000 continue c return end subroutine wrtgry (igray, ngx, ngy, namdmp) c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** ***** c ***** w w rrrr ttttt ggg rrrr y y ***** c ***** w w w r r t g g r r y y ***** c ***** w w w rrrr t g rrrr y y ***** c ***** w w w r r t g gg r r y ***** c ***** w w r r t g g r r y ***** c ***** w w r r t ggg r r y ***** c ***** ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c c we will depend on F77 to automatically allocate/deallocate character*1 ipackd(ngx*ngy) dimension igray (ngx*ngy) character*(*) namdmp c c c c This routine takes the data in the integer*4 array igray and writes c it out to a disk file, in single byte format. c This disk file corresponds one byte per zone to the program data. c c igray = the integer*4 array containing the grayscale integer values. c ngx = the number of zones in a row of the array gray. c ngy = the number of rows in the array gray. c namdmp = an 8-character name for the disk file to be written here. c ipackd = the character*! "packed array. In memory each element of igray c takes 4 bytes, with only the lowest byte used. Hence c the data are NOT contiguous in memory. The data are c contiguous in the array ipackd (which therefore may c be smaller when measured in words.) c c nbytes = ngx*ngy do 3000 i = 1,nbytes ipackd(i) = char(igray(i)) 3000 continue c c c now write it to file namdmp iunit = 34 call ppm98_binopen (namdmp, iunit) call ppm98_binwrite (iunit, ipackd, nbytes) call ppm98_binclose (iunit) c c return end c subroutine runhyd_lr ( xl, yl, ddd, dt, 1 courmx, gamma, smlrho, smalle, smallp, smallu, 2 nx, ny, nbdy) c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** ***** c ***** rrr u u n n h h y y ddd ***** c ***** r r u u nn n h h y y d d ***** c ***** rrr u u n n n hhhhh y d d ***** c ***** rr u u n n n h h y d d ***** c ***** r r u u n nn h h y d d ***** c ***** r r uuu n n h h y ddd ***** c ***** ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension xl( 1-NBDY:NX+NBDY+1) dimension yl( 1-NBDY:NY+NBDY+1) dimension ddd(4,1-nbdy:nx+nbdy,1-nbdy:ny+nbdy) c we are counting on the F77 compiler to automatically allocate/deallocate c c in these arrays, so the dimension nbdy+nx+ny is c a little wasteful, but ensures that enough space has been c allocated. dimension 1 rhonu(1-nbdy:nbdy+nx+ny), ^ pnu(1-nbdy:nbdy+nx+ny), ^ etotnu(1-nbdy:nbdy+nx+ny), ^ u1nu(1-nbdy:nbdy+nx+ny), ^ u2nu(1-nbdy:nbdy+nx+ny) c dimension 1 xxl(1-nbdy:nbdy+nx+ny+1), ^ dddslab(1-nbdy:nbdy+nx+ny,5,6) c c c c nbdy = the number of fake zones needed to implement boundary c conditions so that in the hydrodynamical calculation c all zones can be treated as if they were interior zones. c nx = the number of real zones in the x-dimension for this c tile of the computational domain. c ny = the number of real zones in the y-dimension for this c tile of the computational domain. c iy0x1 = a flag which is 0 in the y-pass and 1 in the x-pass. c c c dt = the time step. c courmx = the maximum value of the Courant number for this tile c of the computational domain. c gamma = the value of the ratio of specific heats, as in the c gamma-law equation of state: p = (gamma-1)*rho*ei c c smlrho = a trivial value of density, used to protect divides. c smalle = a trivial value of energy, used to protect divides. c smallp = a trivial value of pressure, used to protect divides. c smallu = a trivial value of velocity, used to protect divides. c c xl = the x-grid for this tile. c yl = the y-grid for this tile. c rrho = the zone-averaged densities for this tile. c pp = the zone-averaged pressures for this tile. c uux = the zone-averaged x-velocities for this tile. c uuy = the zone-averaged y-velocities for this tile. c ddd = the 3-D array into which the 4 principal 2-D arrays of c basic physical variables are interleaved in order to c optimize memory bandwidth on microprocessor systems. c ddd(1,i,j) = rrho(i,j) c ddd(2,i,j) = eetot(i,j) c ddd(3,i,j) = uux(i,j) c ddd(4,i,j) = uuy(i,j) c c c c 10 continue c courmx = zero maxdim = max (nx, ny) c c c NOW BEGIN LOOP OVER 2 TIME STEPS. c These 2 time steps are symmetrized; they both use the same c time step value, and they perform the x- and y-sweeps in c opposite orders. c iy0x1 = 1 c Set fake zones according to boundary conditions. call bdrys ( xl, yl, ddd, nx, ny, nbdy, iy0x1) PASS(x, y ) 66 format(1h ,i3,1p,4(1x,e10.3)) iy0x1 = 0 c Set fake zones according to boundary conditions. call bdrys ( xl, yl, ddd, nx, ny, nbdy, iy0x1) PASS(y, x ) c iy0x1 = 0 c Set fake zones according to boundary conditions. call bdrys ( xl, yl, ddd, nx, ny, nbdy, iy0x1) PASS(y, x ) c iy0x1 = 1 c Set fake zones according to boundary conditions. call bdrys ( xl, yl, ddd, nx, ny, nbdy, iy0x1) PASS(x, y ) c c c return end c c subroutine bdrys ( xl, yl, ddd, nx, ny, nbdy, iy0x1) c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** ***** c ***** bbbb dddd rrrr y y ssss ***** c ***** b b d d r r y y s s ***** c ***** bbbb d d rrrr y s ***** c ***** b b d d r r y s ***** c ***** b b d d r r y s s ***** c ***** bbbb dddd r r y ssss ***** c ***** ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension xl(-nbdy+1:nx+nbdy+1), yl(-nbdy+1:ny+nbdy+1) dimension ddd(4,1-nbdy:nx+nbdy,1-nbdy:ny+nbdy) c c c c nbdy = the number of fake zones needed to implement boundary c conditions so that in the hydrodynamical calculation c all zones can be treated as if they were interior zones. c nx = the number of real zones in the x-dimension for this c tile of the computational domain. c ny = the number of real zones in the y-dimension for this c tile of the computational domain. c iy0x1 = a flag which is 0 in the y-pass and 1 in the x-pass. c c xl = the x-grid for this tile. c yl = the y-grid for this tile. c c rrho = the zone-averaged densities for this tile. c pp = the zone-averaged pressures for this tile. c uux = the zone-averaged x-velocities for this tile. c uuy = the zone-averaged y-velocities for this tile. c c ddd = the 3-D array into which the 4 principal 2-D arrays of c basic physical variables are interleaved in order to c optimize memory bandwidth on microprocessor systems. c ddd(1,i,j) = rrho(i,j) c ddd(2,i,j) = eet(i,j) c ddd(3,i,j) = uux(i,j) c ddd(4,i,j) = uuy(i,j) c c c c c Apply boundary conditions, Periodic only. c iy1x0 = 1 - iy0x1 nxp1 = nx + 1 ntoset = iy0x1 * nbdy + iy1x0 * 2 c c Periodic boundary conditions: c gridsz = xl(nxp1) - xl(1) do 190 i = 1-nbdy,0 xl(i) = xl(i+nx) - gridsz 190 continue c do 200 i = -ntoset+1,0 do 200 k = 1,ny ddd(1,i,k) = ddd(1,i+nx,k) ddd(2,i,k) = ddd(2,i+nx,k) ddd(3,i,k) = ddd(3,i+nx,k) ddd(4,i,k) = ddd(4,i+nx,k) 200 continue c do 290 i = nxp1+1,nxp1+nbdy xl(i) = xl(i-nx) + gridsz 290 continue c do 300 i = nxp1,nx+ntoset do 300 k = 1,ny ddd(1,i,k) = ddd(1,i-nx,k) ddd(2,i,k) = ddd(2,i-nx,k) ddd(3,i,k) = ddd(3,i-nx,k) ddd(4,i,k) = ddd(4,i-nx,k) 300 continue c 900 continue c c nyp1 = ny + 1 ntoset = iy1x0 * nbdy + iy0x1 * 2 c c Periodic boundary conditions: c gridsz = yl(nyp1) - yl(1) do 1190 i = 1-nbdy,0 yl(i) = yl(i+ny) - gridsz 1190 continue c do 1200 i = -ntoset+1,0 do 1200 k = -nbdy+1,nx+nbdy ddd(1,k,i) = ddd(1,k,i+ny) ddd(2,k,i) = ddd(2,k,i+ny) ddd(3,k,i) = ddd(3,k,i+ny) ddd(4,k,i) = ddd(4,k,i+ny) 1200 continue c do 1290 i = nyp1+1,nyp1+nbdy yl(i) = yl(i-ny) + gridsz 1290 continue c do 1300 i = nyp1,ny+ntoset do 1300 k = -nbdy+1,nx+nbdy ddd(1,k,i) = ddd(1,k,i-ny) ddd(2,k,i) = ddd(2,k,i-ny) ddd(3,k,i) = ddd(3,k,i-ny) ddd(4,k,i) = ddd(4,k,i-ny) 1300 continue c c c return end subroutine eos( rho, e_tot, u1, u2, p, ce, gamma, n ) c c c This routine applies the equation of state, assuin a c gamma law EOS: p(i) = ( gamma - 1.0 ) * rho(i) * ei(I) c ei(i) = e_tot(i) - 0.5*( u1(i)**2 + u2(i)**2) c c where c rho(i) is the zone averaged density, c p(i) is the zone averaged pressure, c ei(i) is the zone averaged total energy, c ei(i) is the zone averaged internal energy, c u1(i) is the zone averaged velocity-1, c u2(i) is the zone averaged velocity-2, c ce(i) is the zone averaged Eulerian sound speed c c gamma is the ration of specific heats of the gas. c c n is the number of zones c c c ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension rho(n), e_tot(n), u1(n), u2(n), p(n), ce(n) gamma1 = gamma - one ggamma1 = gamma * gamma1 do i=1,n ek = 0.5 * ( u1(i)**2 + u2(i)**2 ) ei = e_tot(i) - ek ei = max ( zero, ei ) p(i) = gamma1 * rho(i) * ei ce(i) = sqrt( ggamma1 * ei ) enddo return end subroutine wrtimg4 (t1, t2, t3, t4, acharc, acharc2, nx, ny, iimage ) c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** ***** c ***** w w rrrr ttttt iii m m ggg ***** c ***** w w w r r t i m m g g ***** c ***** ww ww rrrr t i mm mm g ***** c ***** ww ww r r t i mm mm g gg ***** c ***** w w r r t i m m m g g ***** c ***** w w r r t iii m m m ggg ***** c ***** ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension t1(nx,ny) dimension t2(nx,ny) dimension t3(nx,ny) dimension t4(nx,ny) parameter (nngx=IQX) parameter (nngy=IQY) parameter (ng=nngx*nngy) dimension igray(ng) dimension ibiggray(4*ng) character*8 namdmp character*1 acharc character*1 acharc2 c c t1 = a 2-D temparay array containing the data to write to disk c c final combined images will be: c ---------- c |t1 | t2| c | | | c ---------- c | | | c |t3 | t4| c ---------- c nx = the x dimension of the array t. c ny = the y dimension of the array t. c iimage = the sequential number of the image 0->whatever c igray = the 2-D array of grayscale values (ranging from 0 to 255) c ibiggray = the 2-D combined array of 4 images c acharc = the first character in the names of all images files c acharc2 = the second character in the names of all images files c c all binary image files are named 'abim####' where a = acharc, c b = acharc2, and #### is the sequential number 0000 to 9999 c namdmp = the 8-character name of the image file to be created. c c c Construct name of image file and store in NAMDMP: c write (namdmp,1004) acharc, acharc2, iimage 1004 format (a1,a1,'im',i4.4) c c c Scale array T so that it ranges from 0 to 255, storing the result in c the integer array IGRAY. Also returns the minimum (TMPMIN) and maximum c (TMPMAX) of array T. c call makgry (t1, igray, tmpmin, tmpmax, nx, ny) do j=1,ny do i=1,nx ii = i + nx*(j-1) iii = i + 2*nx*(j-1) ibiggray(iii) = igray(ii) enddo enddo call makgry (t2, igray, tmpmin, tmpmax, nx, ny) do j=1,ny do i=1,nx ii = i + nx*(j-1) iii = nx+i + 2*nx*(j-1) ibiggray(iii) = igray(ii) enddo enddo call makgry (t3, igray, tmpmin, tmpmax, nx, ny) do j=1,ny do i=1,nx ii = i + nx*(j-1) iii = i + 2*nx*(j+ny-1) ibiggray(iii) = igray(ii) enddo enddo call makgry (t4, igray, tmpmin, tmpmax, nx, ny) do j=1,ny do i=1,nx ii = i + nx*(j-1) iii = nx+i + 2*nx*(j+ny-1) ibiggray(iii) = igray(ii) enddo enddo c c Now write it to disk. c call wrtgry (ibiggray, 2*nx, 2*ny, namdmp) c return end