define(IQX,128) define(IQY,128) define(IQZ,128) define(IQ,128) define(QNBDY,5) define(PREFIX,L) 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,z ) (without the "_") c in which case, all instances of "$1" will be replaced with "x", and c all "$2" by "y", and all "$3" by "z". c c Why has this been done? Believe it or not, as an attempt at simplification. c In order to update any given row, PPM requires information from up to c two rows above below, fore, and aft of the row under intense scrutiny. c c Consider first just above and below. c c This information must be from the same computational time, however, c as the computation progresses, the large storage array 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 five rows deep (a slab) will c be created. Consider an X-pass, then the slab would be dimensioned c at least (1-nbdy:nx+nbdy,,-1,ny+2,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 planes (islab1bb, islab1b, islab1, islab1t, c islab1tt) c (for 2 planes aft, 1 plane aft, the row, one row fore, 2 rows fore) c which would have the values (1,2,3,4,5). After the call to do_ppm and the c updated values stored in DDD, the counters are cylcled to (2,3,4,5,1) c and the relavent required data from DDD is written into the row pointed to c by the new value of islabtt. 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, u3 are the y- and c z-velocitiets, resp. while during a Y-pass, u1 is the y-velocity and u2, u3 c are the z- and x-velocities, resp. One should maintain a right-hand c coordinate system. I don't believe that this is required (now) but c who knows what future perils await thee. c c c This reduces both the scratch space requirement and the amount of c data copying required. Believe it or not. c c the xvelocity is element #3 in DDD c the yvelocity is element #4 in DDD c the zvelocity is element #5 in DDD kvx = 3 kvy = 4 kvz = 5 do i = 1-nbdy, n$1+nbdy+1 xxl(i) = $1l(i) enddo c c initialize the first slab c c islab = 5 is done below at plane islabtt c c ((do 100) do islab = 1, 4 i$3 = islab-2 do i$2 = 1-nbdy,n$2+nbdy do i$1 = 1-nbdy, n$1+nbdy dddslab(i$1, i$2,islab,1) = ddd( 1, ix, iy,iz) dddslab(i$1, i$2,islab,2) = ddd( 2, ix, iy,iz) dddslab(i$1, i$2,islab,3) = ddd(kv$1, ix, iy,iz) dddslab(i$1, i$2,islab,4) = ddd(kv$2, ix, iy,iz) dddslab(i$1, i$2,islab,5) = ddd(kv$3, ix, iy,iz) enddo enddo enddo c (do 100)) islabbb = 0 islabb = 1 islab = 2 islabt = 3 islabtt = 4 c c ((do 410) do i$3 = 1, n$3 c (slabbb) two in front (slabtt), etc. c islabbb = 1+mod(i$3-1,5) c islabb = 1+mod(i$3, 5) c islab = 1+mod(i$3+1,5) c islabt = 1+mod(i$3+2,5) c islabtt = 1+mod(i$3+3,5) islabbb = 1+mod(islabbb,5) islabb = 1+mod(islabb, 5) islab = 1+mod(islab, 5) islabt = 1+mod(islabt, 5) islabtt = 1+mod(islabtt,5) c At this point, we are ALWAYS able to fill plane islabtt from ddd l$3 = i$3+2 do l$2 = 1-nbdy, n$2+nbdy do l$1 = 1-nbdy, n$1+nbdy dddslab(l$1, l$2,islabtt,1) = ddd( 1, lx, ly,lz) dddslab(l$1, l$2,islabtt,2) = ddd( 2, lx, ly,lz) dddslab(l$1, l$2,islabtt,3) = ddd(kv$1, lx, ly,lz) dddslab(l$1, l$2,islabtt,4) = ddd(kv$2, lx, ly,lz) dddslab(l$1, l$2,islabtt,5) = ddd(kv$3, lx, ly,lz) enddo enddo c ((do 400) do i$2 = 1, n$2 c ip = 1-nbdy $1courmx = zero call do_ppm_de0_3dc_gamma( xl(ip), ^ dddslab(ip,i$2-2,islab,1), ^ dddslab(ip,i$2,islab,1), ^ dddslab(ip,i$2+2,islab,1), ^ dddslab(ip,i$2,islabbb,1), ^ dddslab(ip,i$2,islabtt,1), ^ dddslab(ip,i$2-2,islab,2), ^ dddslab(ip,i$2,islab,2), ^ dddslab(ip,i$2+2,islab,2), ^ dddslab(ip,i$2,islabbb,2), ^ dddslab(ip,i$2,islabtt,2), ^ dddslab(ip,i$2-2,islab,3), ^ dddslab(ip,i$2,islab,3), ^ dddslab(ip,i$2+2,islab,3), ^ dddslab(ip,i$2,islabbb,3), ^ dddslab(ip,i$2,islabtt,3), ^ dddslab(ip,i$2-2,islab,4), ^ dddslab(ip,i$2-1,islab,4), ^ dddslab(ip,i$2,islab,4), ^ dddslab(ip,i$2+1,islab,4), ^ dddslab(ip,i$2+2,islab,4), ^ dddslab(ip,i$2,islabbb,4), ^ dddslab(ip,i$2,islabtt,4), ^ dddslab(ip,i$2,islabbb,5), ^ dddslab(ip,i$2,islabb,5), ^ dddslab(ip,i$2,islab,5), ^ dddslab(ip,i$2,islabt,5), ^ dddslab(ip,i$2,islabtt,5), ^ dddslab(ip,i$2-2,islab,5), ^ dddslab(ip,i$2+2,islab,5), ^ rhonu(ip), pnu(ip), u1nu(ip), u2nu(ip), u3nu(ip), ^ gamma, dt, smlro, smallp, smallu, smalle, ^ $1courmx, n$1, nbdy ) c c print *,' i$2, i$3, $1courmx=', i$2, i$3, $1courmx c courmx = max ( courmx, $1courmx ) c do i$1 = 1, n$1 ddd( 1,ix,iy,iz) = rhonu(i$1) ddd( 2,ix,iy,iz) = pnu(i$1) ddd(kv$1,ix,iy,iz) = u1nu(i$1) ddd(kv$2,ix,iy,iz) = u2nu(i$1) ddd(kv$3,ix,iy,iz) = u3nu(i$1) enddo c enddo c (do 400)) enddo c (do 410)) c c c c ) c) C############################################################### c======================================================================= c= = c= It is assumed that this file will be preprocessed by m4 = c= = c======================================================================= c c single precision: c m4 -B6666 -Difdouble=0 Driverlagremap3dg0c.EXA > driver.f c c double precision: c m4 -B6666 -Difdouble=1 Driverlagremap3dg0c.EXA > driver.f c c c c======================================================================= 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 IQZ gives the number of zones of a physical column c (z-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, IQZ 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 Driverdirect3dg0c * 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 direct Eulerian style c using routines from the PPMLIB library. Directional c splitting is used to approximate the 3-dimensional c Euler equations by applying a 1-dimensional operator c to treat gradients in the x-, y- and z- directions c separately. The inital problem consists of a three- c dimensional shock tube (discontinuities in x-, y- and c z-directions dividing the initial state into eqight c states 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 3 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 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 5 "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 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 Four images are included in each file in a 2x2 matrix. c These images are density, pressure and x- and y- velocities. c Thus each file is 2*NX x 2*NY bytes giving a one pixel per c zone ratio when displayed with the XRAZ routines c available from LCSE and elsewhere. c c C program Driverdirect2dg0c c c parameter (NXZONES = IQX) parameter (NYZONES = IQY) parameter (NZZONES = IQZ) parameter (MAXSTEPS = 500) parameter (EOSGAM = 1.4) parameter (DT = 0.00200) parameter (XMAX = 1.0) parameter (YMAX = 1.0) parameter (ZMAX = 1.0) parameter (N2DUMP = 4) 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 NZZONES = number of real computational zones in the z-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 ZMAX = the righthand edge y-value of the grid. c The z-grid runs from (0,ZMAX) 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 = QNBDY) 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 zl( 1-NBDY:NYZONES+NBDY+1) dimension ddd( 1:5, ^ 1-NBDY:NXZONES+NBDY, ^ 1-NBDY:NYZONES+NBDY, ^ 1-NBDY:NZZONES+NBDY ) c c c c ddd = the 4-D array into which the 5 principal 3-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 zl = the locations of the front zone interfaces c ddd = the 4-D array into which the 5 principal 3-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) = pp(i,j) c ddd(3,i,j) = uux(i,j) c ddd(4,i,j) = uuy(i,j) c ddd(5,i,j) = uuz(i,j) c rrho = the zone-averaged densities c pp = the zone-averaged pressures c uux = the zone-averaged x-velocities c uuy = the zone-averaged y-velocities c uuz = the zone-averaged z-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 zlimit = ZMAX nx = nxzones ny = nyzones nz = nzzones c call setgrid( xl, yl, zl, xlimit, ylimit, zlimit, ^ nx, ny, nz, nbdy ) c c c initialize the hydrodynamic quantities: c call set3dshock( xl, yl, zl, ddd, 5, nx, ny, nz, 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, zl, ddd, gamma, dtime, smlrho, smalle, smallp, ^ smallu, numdmp, ndump, tdump, nstop, tstop, nx, ny, ^ nz, nbdy) stop end subroutine cntrol ( xl, yl, zl, ddd, gamma, dtime, smlrho, smalle, ^ smallp, smallu, numdmp, ndump, tdump, nstop, tstop, nx, ny, nz, 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 zl( 1-NBDY:NZ+NBDY+1) dimension ddd( 1:5, 1-NBDY:NX+NBDY, 1-NBDY:NY+NBDY, 1-NBDY:NZ+NBDY ) c c c This routine runs the entire hydro calculation. c c c c ddd = the 4-D array into which the 5 principal 3-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 nz = the number of real zones in the z-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 zl = the locations of the front zone interfaces c ddd = the 3-D array into which the 5 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) = pp(i,j) c ddd(3,i,j) = uux(i,j) c ddd(4,i,j) = uuy(i,j) c ddd(5,i,j) = uuz(i,j) 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 uuz = the zone-averaged z-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 dump3d (ddd, numdmp, nx, ny, nz, 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 runhyd 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 RUNHYD3 performs a symmetrized pair of passes x-y-z-z-y-x which c equals two complete timesteps. c call runhyd3 ( xl, yl, zl, ddd, dtime, courmx, gamma, 1 smlrho, smalle, smallp, smallu, numdmp, 2 nx, ny, nz, 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 dump3d (ddd, numdmp, nx, ny, nz, 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, zl, xmax, ymax, zmax, nx, ny, nz, 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 zl = the z-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 zmax = The real (non boundary) portion of the grid is c ZZL(1)=0.0 to ZZL(NZ+1)=ZMAX. c c nx = number of real zones in in the x-direction c ny = number of real zones in in the y-direction c nz = number of real zones in in the z-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) dimension zl(1-nbdy:nz+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 deez = zmax / float(NZ) DO I=1,NZ+1 zl(i) = deez*float(i-1) ENDDO return end c subroutine runhyd3 ( xl, yl, zl, ddd, dt, 1 courmx, gamma, smlrho, smalle, smallp, smallu, 2 numdmp, nx, ny, nz, 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 zl( 1-NBDY:NZ+NBDY+1) dimension ddd(5,1-nbdy:nx+nbdy,1-nbdy:ny+nbdy,1-nbdy:nz+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-QNBDY:QNBDY+IQ), ^ pnu(1-QNBDY:QNBDY+IQ), ^ u1nu(1-QNBDY:QNBDY+IQ), ^ u2nu(1-QNBDY:QNBDY+IQ), ^ u3nu(1-QNBDY:QNBDY+IQ), ^ crmx(1-QNBDY:QNBDY+IQ) c dimension 1 xxl(1-nbdy:nbdy+nx+ny+1), ^ dddslab(1-QNBDY:QNBDY+IQ,1-QNBDY:QNBDY+IQ,5,5) 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 nz = the number of real zones in the z-dimension for this c tile of the computational domain. 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 zl = the z-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 uuz = the zone-averaged y-velocities for this tile. c ddd = the 3-D array into which the 5 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) = pp(i,j) c ddd(3,i,j) = uux(i,j) c ddd(4,i,j) = uuy(i,j) c ddd(5,i,j) = uuz(i,j) c c c c 10 continue c courmx = zero maxdim = max (nx, ny,nz) 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 sweeps in c opposite orders. c c Set fake zones according to boundary conditions. call bdrys ( xl, yl, zl, ddd, 5, 5, 5, nx, ny, nz, nbdy ) PASS(x, y, z ) print *,' pass x, y, z 1 courmx = ',courmx 66 format(1h ,i3,1p,4(1x,e10.3)) c Set fake zones according to boundary conditions. call bdrys ( xl, yl, zl, ddd, 5, 5, 5, nx, ny, nz, nbdy ) PASS(y, z, x ) print *,' pass y, z, x 1 courmx = ',courmx c c Set fake zones according to boundary conditions. call bdrys ( xl, yl, zl, ddd, 5, 5, 5, nx, ny, nz, nbdy ) PASS(z, x, y ) print *,' pass z, x, y 1 courmx = ',courmx c c Set fake zones according to boundary conditions. call bdrys ( xl, yl, zl, ddd, 5, 5, 5, nx, ny, nz, nbdy ) PASS(z, x, y ) print *,' pass z, x, y 2 courmx = ',courmx c c Set fake zones according to boundary conditions. call bdrys ( xl, yl, zl, ddd, 5, 5, 5, nx, ny, nz, nbdy ) PASS(y, z, x ) print *,' pass y, z, x 2 courmx = ',courmx c c Set fake zones according to boundary conditions. call bdrys ( xl, yl, zl, ddd, 5, 5, 5, nx, ny, nz, nbdy ) PASS(x, y, z ) print *,' pass x, y, z 2 courmx = ',courmx c return end c c subroutine bdrys ( xl, yl, zl, ddd, nxtoset, nytoset, nztoset, ^ nx, ny, nz, nbdy ) 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 Zl(-nbdy+1:nZ+nbdy+1) dimension ddd(5,1-nbdy:nx+nbdy, ^ 1-nbdy:ny+nbdy, ^ 1-nbdy:nz+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 nxtoset= the number of boundary zones to set at each end in c the x-dimension for this pass. c nytoset= the number of boundary zones to set at each end in c the y-dimension for this pass. c nztoset= the number of boundary zones to set at each end in c the z-dimension for this pass. 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 nz = the number of real zones in the z-dimension for this c tile of the computational domain. c c xl = the x-grid for this tile. c yl = the y-grid for this tile. c zl = the z-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 uuz = the zone-averaged z-velocities for this tile. c c ddd = the 4-D array into which the 5 principal 3-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) = pp(i,j) c ddd(3,i,j) = uux(i,j) c ddd(4,i,j) = uuy(i,j) c ddd(5,i,j) = uuz(i,j) c c c c c Apply boundary conditions, reflecting only. c c c c Reflecting boundary conditions: c do i = 1-nxtoset,0 xl(i) = 2.0*xl(1) - xl(2-i) enddo do i = nx+2,nx+nxtoset+1 xl(i) = 2.0*xl(nx+1) - xl(2*nx+2-i) enddo c do i = 1-nytoset,0 yl(i) = 2.0*yl(1) - yl(2-i) enddo do i = ny+2,ny+nytoset+1 yl(i) = 2.0*yl(ny+1) - yl(2*ny+2-i) enddo c do i = 1-nztoset,0 zl(i) = 2.0*zl(1) - zl(2-i) enddo do i = nz+2,nz+nztoset+1 zl(i) = 2.0*zl(nz+1) - zl(2*nz+2-i) enddo c c c reflecting boundary conditions: c c do 1200 k = 1-nztoset,0 do 1200 j = 1,ny do 1200 i = 1,nx ddd(1,i,j,k) = ddd(1,i,j,1-k) ddd(2,i,j,k) = ddd(2,i,j,1-k) ddd(3,i,j,k) = ddd(3,i,j,1-k) ddd(4,i,j,k) = ddd(4,i,j,1-k) ddd(5,i,j,k) = -ddd(5,i,j,1-k) 1200 continue c do 1210 k = nz+1,nz+nztoset do 1210 j = 1,ny do 1210 i = 1,nx ddd(1,i,j,k) = ddd(1,i,j,2*nz+1-k) ddd(2,i,j,k) = ddd(2,i,j,2*nz+1-k) ddd(3,i,j,k) = ddd(3,i,j,2*nz+1-k) ddd(4,i,j,k) = ddd(4,i,j,2*nz+1-k) ddd(5,i,j,k) =-ddd(5,i,j,2*nz+1-k) 1210 continue c do 1300 k = 1-nztoset,nz+nztoset do 1300 j = 1-nytoset,0 do 1300 i = 1,nx ddd(1,i,j,k) = ddd(1,i,1-j,k) ddd(2,i,j,k) = ddd(2,i,1-j,k) ddd(3,i,j,k) = ddd(3,i,1-j,k) ddd(4,i,j,k) =-ddd(4,i,1-j,k) ddd(5,i,j,k) = ddd(5,i,1-j,k) 1300 continue c do 1310 k = 1-nztoset,nz+nztoset do 1310 j = ny+1,ny+nytoset do 1310 i = 1,nx ddd(1,i,j,k) = ddd(1,i,2*ny+1-j,k) ddd(2,i,j,k) = ddd(2,i,2*ny+1-j,k) ddd(3,i,j,k) = ddd(3,i,2*ny+1-j,k) ddd(4,i,j,k) =-ddd(4,i,2*ny+1-j,k) ddd(5,i,j,k) = ddd(5,i,2*ny+1-j,k) 1310 continue c do 1400 k = 1-nztoset,nz+nztoset do 1400 j = 1-nytoset,ny+nytoset do 1400 i = 1-nxtoset,0 ddd(1,i,j,k) = ddd(1,1-i,j,k) ddd(2,i,j,k) = ddd(2,1-i,j,k) ddd(3,i,j,k) =-ddd(3,1-i,j,k) ddd(4,i,j,k) = ddd(4,1-i,j,k) ddd(5,i,j,k) = ddd(5,1-i,j,k) 1400 continue c do 1410 k = 1-nztoset,nz+nztoset do 1410 j = 1-nytoset,ny+nytoset do 1410 i = nx+1,nx+nxtoset ddd(1,i,j,k) = ddd(1,2*nx+1-i,j,k) ddd(2,i,j,k) = ddd(2,2*nx+1-i,j,k) ddd(3,i,j,k) =-ddd(3,2*nx+1-i,j,k) ddd(4,i,j,k) = ddd(4,2*nx+1-i,j,k) ddd(5,i,j,k) = ddd(5,2*nx+1-i,j,k) 1410 continue c c c c return end subroutine dump3d (ddd, numdmp, nx, ny, nz, nbdy ) c c This will create 5 binary image files scaled to one byte per pixel c (0-155), one pixel per zone, arranges as a "brick of bytes" suitable c for views by "bob." c c The naming convention is PREFIX num2.2 'm' num4.4 c c where PREFIX is a one character alpha numeric identifier, c num2.2 is a tow digit identifier: c 00 = density c 01 = pressure c 02 = x-velocity c 03 = y-velocity c 04 = z-velocity c num4.4 is the 4 digit sequence number (0000 - 9999) 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 character*1 igray(nx,ny) character*8 namdmp c c c ddd = the 4-D array into which the 5 principal 3-D arrays of c basic physical variables are interleaved in order to c optimize memory bandwidth on microprocessor systems. c ddd(1,i,j,k) = rrho(i,j,k) c ddd(2,i,j,k) = pp(i,j,k) c ddd(3,i,j,k) = uux(i,j,k) c ddd(4,i,j,k) = uuy(i,j,k) c ddd(5,i,j,k) = uuz(i,j,k) c rrho = the zone-averaged densities. c pp = the zone-averaged pressures. c uux = the zone-averaged x-velocities. c uuy = the zone-averaged y-velocities. c uuz = the zone-averaged z-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 nz = number of real computational zones in the z-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:5, 1-NBDY:NX+NBDY, 1-NBDY:NY+NBDY, 1-NBDY:NZ+NBDY ) c c acharc = 'PREFIX' c c determine individual min/max for scaling c rhomin = ddd(1,1,1,1) rhomax = ddd(1,1,1,1) pmin = ddd(2,1,1,1) pmax = ddd(2,1,1,1) uxmin = ddd(3,1,1,1) uxmax = ddd(3,1,1,1) uymin = ddd(4,1,1,1) uymax = ddd(4,1,1,1) uzmin = ddd(5,1,1,1) uzmax = ddd(5,1,1,1) do 1000 k = 1,nz do 1000 j = 1,ny do 1000 i = 1,nx ii = ii + 1 rhomin = min( rhomin,ddd(1,i,j,k)) rhomax = max( rhomax,ddd(1,i,j,k)) pmin = min( pmin,ddd(2,i,j,k)) pmax = max( pmax,ddd(2,i,j,k)) uxmin = min( uxmin,ddd(3,i,j,k)) uxmax = max( uxmax,ddd(3,i,j,k)) uymin = min( uymin,ddd(4,i,j,k)) uymax = max( uymax,ddd(4,i,j,k)) uzmin = min( uzmin,ddd(5,i,j,k)) uzmax = max( uzmax,ddd(5,i,j,k)) 1000 continue c c create scaling factors c denom = (rhomax - rhomin) if( denom .ne. 0.0 ) denom = 1.0/denom rhofac = 255.0 * denom denom = (pmax - pmin) if( denom .ne. 0.0 ) denom = 1.0/denom pfac = 255.0 * denom denom = (uxmax - uxmin) if( denom .ne. 0.0 ) denom = 1.0/denom uxfac = 255.0 * denom denom = (uymax - uymin) if( denom .ne. 0.0 ) denom = 1.0/denom uyfac = 255.0 * denom denom = (uzmax - uzmin) if( denom .ne. 0.0 ) denom = 1.0/denom uzfac = 255.0 * denom c c Construct name of image file and store in NAMDMP: write (namdmp,1004) acharc, numdmp 1004 format (a1,'01m',i4.4) c now write it to file namdmp iunit = 34 nbytes = nx*ny call ppm98_binopen (namdmp, iunit) do k = 1,nz do j = 1,ny do i = 1,nx rthyng = rhofac * (ddd(1,i,j,k) - rhomin) igray(i,j) = char(int(rthyng)) enddo enddo call ppm98_binwrite (iunit, igray, nbytes) enddo call ppm98_binclose (iunit) c c c Construct name of image file and store in NAMDMP: write (namdmp,1005) acharc, numdmp 1005 format (a1,'02m',i4.4) c now write it to file namdmp iunit = 34 call ppm98_binopen (namdmp, iunit) do k = 1,nz do j = 1,ny do i = 1,nx rthyng = pfac * (ddd(2,i,j,k) - pmin) igray(i,j) = char(int(rthyng)) enddo enddo call ppm98_binwrite (iunit, igray, nbytes) enddo call ppm98_binclose (iunit) c c c Construct name of image file and store in NAMDMP: write (namdmp,1006) acharc, numdmp 1006 format (a1,'03m',i4.4) c now write it to file namdmp iunit = 34 call ppm98_binopen (namdmp, iunit) do k = 1,nz do j = 1,ny do i = 1,nx rthyng = uxfac * (ddd(3,i,j,k) - uxmin) igray(i,j) = char(int(rthyng)) enddo enddo call ppm98_binwrite (iunit, igray, nbytes) enddo call ppm98_binclose (iunit) c c c Construct name of image file and store in NAMDMP: write (namdmp,1007) acharc, numdmp 1007 format (a1,'04m',i4.4) c now write it to file namdmp iunit = 34 call ppm98_binopen (namdmp, iunit) do k = 1,nz do j = 1,ny do i = 1,nx rthyng = uyfac * (ddd(4,i,j,k) - uymin) igray(i,j) = char(int(rthyng)) enddo enddo call ppm98_binwrite (iunit, igray, nbytes) enddo call ppm98_binclose (iunit) c c c Construct name of image file and store in NAMDMP: write (namdmp,1008) acharc, numdmp 1008 format (a1,'05m',i4.4) c now write it to file namdmp iunit = 34 call ppm98_binopen (namdmp, iunit) do k = 1,nz do j = 1,ny do i = 1,nx rthyng = uzfac * (ddd(5,i,j,k) - uzmin) igray(i,j) = char(int(rthyng)) enddo enddo call ppm98_binwrite (iunit, igray, nbytes) enddo call ppm98_binclose (iunit) c c numdmp = numdmp + 1 c c return end subroutine set3dshock( xl, yl, zl, ddd, nvar, nx, ny, nz, nbdy) c ************************************************************************** c * THIS CODE IS D36PPM3D * c * THIS CODE WAS WRITTEN BY: B. Kevin Edgar ALL RIGHTS RESERVED * c ************************************************************************** dimension xl(1-nbdy:nx+nbdy+1) dimension yl(1-nbdy:ny+nbdy+1) dimension zl(1-nbdy:nz+nbdy+1) dimension ddd(1:nvar,1-nbdy:nx+nbdy, ^ 1-nbdy:ny+nbdy, ^ 1-nbdy:nz+nbdy) c do k=1-nbdy,nz+nbdy do j=1-nbdy,ny+nbdy do i=1-nbdy,nx+nbdy do l=1,nvar ddd(l,i,j,k) = 1.0 enddo enddo enddo enddo xdisc = 0.49 ydisc = 0.49 zdisc = 0.49 xlength = xl(nx+1) - xl(1) ylength = yl(ny+1) - yl(1) zlength = zl(nz+1) - zl(1) xdisc = xl(1) + xdisc * xlength ydisc = yl(1) + ydisc * ylength zdisc = zl(1) + zdisc * zlength c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c This routine sets the pressure, density, and velocity fields for c c a BKE 3d shocktube. c c----------------------------------------------------------------------c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccc rhoo8 = 2.5 po8 = 2.5 vxo8 = 0.0 vyo8 = 0.0 vzo8 = 0.0 rhoo7 = 1.25 po7 = 1.25 vxo7 = 0.0 vyo7 = 0.0 vzo7 = 0.0 rhoo6 = 0.125 po6 = 0.125 vxo6 = 0.0 vyo6 = 0.0 vzo6 = 0.0 rhoo5 = 0.635 po5 = 0.635 vxo5 = 0.0 vyo5 = 0.0 vzo5 = 0.0 c rhoo1 = 0.375 po1 = 0.375 vxo1 = 0.0 vyo1 = 0.0 vzo1 = 0.0 rhoo2 = 0.075 po2 = 0.075 vxo2 = 0.0 vyo2 = 0.0 vzo2 = 0.0 rhoo3 = 0.75 po3 = 0.75 vxo3 = 0.0 vyo3 = 0.0 vzo3 = 0.0 rhoo4 = 1.5 po4 = 1.5 vxo4 = 0.0 vyo4 = 0.0 vzo4 = 0.0 c c y c ^ c | c ******************* c * 6 * 5 ** c * * * * c ******************* * c * * ** * c * * * * * c ******************* * ** c * * * ** * c * 2 * 1 * * * c * * * **8 * c * * ** * *-->x c ******************* * * c * * * ** c * 3 * 4 * * c * * * * c * * ** c ******************* c / c z c do 611 k=1, nz if ( zl(k) .ge. zdisc ) then rrhoq1 = rhoo1 rrhoq2 = rhoo2 rrhoq3 = rhoo3 rrhoq4 = rhoo4 ppq1 = po1 ppq2 = po2 ppq3 = po3 ppq4 = po4 vvxq1 = vxo1 vvxq2 = vxo2 vvxq3 = vxo3 vvxq4 = vxo4 vvyq1 = vyo1 vvyq2 = vyo2 vvyq3 = vyo3 vvyq4 = vyo4 vvzq1 = vzo1 vvzq2 = vzo2 vvzq3 = vzo3 vvzq4 = vzo4 else rrhoq1 = rhoo5 rrhoq2 = rhoo6 rrhoq3 = rhoo7 rrhoq4 = rhoo8 ppq1 = po5 ppq2 = po6 ppq3 = po7 ppq4 = po8 vvxq1 = vxo5 vvxq2 = vxo6 vvxq3 = vxo7 vvxq4 = vxo8 vvyq1 = vyo5 vvyq2 = vyo6 vvyq3 = vyo7 vvyq4 = vyo8 vvzq1 = vzo5 vvzq2 = vzo6 vvzq3 = vzo7 vvzq4 = vzo8 endif do 611 j=1, ny if ( yl(j) .ge. ydisc ) then rhol = rrhoq2 pl = ppq2 vxl = vvxq2 vyl = vvyq2 vzl = vvzq2 rhor = rrhoq1 pr = ppq1 vxr = vvxq1 vyr = vvyq1 vzr = vvzq1 else rhol = rrhoq3 pl = ppq3 vxl = vvxq3 vyl = vvyq3 vzl = vvzq3 rhor = rrhoq4 pr = ppq4 vxr = vvxq4 vyr = vvyq4 vzr = vvzq4 endif do 611 i=1, nx if ( xl(i) .le. xdisc ) then ddd(1,i,j,k) = rhor ddd(2,i,j,k) = pr ddd(3,i,j,k) = vxr ddd(4,i,j,k) = vyr ddd(5,i,j,k) = vzr else ddd(1,i,j,k) = rhol ddd(2,i,j,k) = pl ddd(3,i,j,k) = vxl ddd(4,i,j,k) = vyl ddd(5,i,j,k) = vzl endif 611 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c return end