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 Driverdirect2dg0c_v01.EXA > driver.f c c double precision: c m4 -B6666 -Difdouble=1 Driverdirect2dg0c_v01.EXA > driver.f c c 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 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 row 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 call ifelse(ifdouble,1,d_, ) do_ppmde0_2dc_gamma( xxl(ip), ^ dddslab(ip,islabbb,1), dddslab(ip,islab,1), dddslab(ip,islabtt,1), ^ dddslab(ip,islabbb,2), dddslab(ip,islab,2), dddslab(ip,islabtt,2), ^ dddslab(ip,islabbb,3), dddslab(ip,islab,3), dddslab(ip,islabtt,3), ^ dddslab(ip,islabbb,4), dddslab(ip,islabb,4), dddslab(ip,islab,4), ^ dddslab(ip,islabt,4), dddslab(ip,islabtt,4), ^ rhonu(ip), pnu(ip), u1nu(ip), u2nu(ip), ^ gamma, dt, smlrho, smallp, smallu, ^ smalle, $1courmx, n$1, nbdy) c call do_ppmde0_2dc_gamma ( xl(L) c ^ rrho(L,k-2), rrho(L,k), rrho(L,k+2), c ^ pp(L,k-2), pp(L,k), pp(L,k+2), c ^ uux(L,k-2), uux(L,k), uux(L,k+2), c ^ uuy(L,k-2), uuy(L,k-1), uuy(L,k), uuy(L,k+1), uuy(L,k+2), c ^ rhonu(L), pnu(L), uxnu(L), uynu(L), c ^ gamma, dt, smlrho, smallp, smallu, c ^ smalle, xcourmx, nx, nbdy) c call ifelse(ifdouble,1,d_, ) do_ppmlr0_2dc_gamma ( xxl(ip), c ^ dddslab(ip,islabbb,1), dddslab(ip,islab,1), dddslab(ip,islabtt,1), c ^ dddslab(ip,islabbb,2), dddslab(ip,islab,2), dddslab(ip,islabtt,2), c ^ dddslab(ip,islabbb,3), dddslab(ip,islab,3), dddslab(ip,islabtt,3), c ^ dddslab(ip,islabbb,4), dddslab(ip,islabb,4), dddslab(ip,islab,4), c ^ dddslab(ip,islabt,4), dddslab(ip,islabtt,4), c ^ rhonu(ip), pnu(ip), u1nu(ip), u2nu(ip), c ^ gamma, dt, smlrho, smallp, smallu, c ^ 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) = pnu(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= It is assumed that this file will be preprocessed by m4 = c= = c======================================================================= c c single precision: c m4 -B6666 -Difdouble=0 Driverdirect2dg0c_v01.EXA > driver.f c c double precision: c m4 -B6666 -Difdouble=1 Driverdirect2dg0c_v01.EXA > driver.f c c c c======================================================================= define(IQX,256) define(IQY,256) define(IQ,256) define(PREFIX,D) 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 Driverdirect2dg0c_v01 * c * It is derived from Driverdirect2dg0c_v00 * 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 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 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_v01 c c parameter (NXZONES = IQX) parameter (NYZONES = IQY) parameter (EOSGAM = 1.4) parameter (DT = 0.001100) 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 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 = 5) 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 ) character*1 achar c c c 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 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) = pp(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 pp = the zone-averaged pressures 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 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 initialize the grid: xlimit = XMAX ylimit = YMAX c call setgrid( xl, yl, xlimit, ylimit, nxzones, nyzones, nbdy ) c print *,' Enter 1 if you intend to restart from file ''RESTART''' print *,' Enter anything else to start the problem with the initial' print *,' conditions.' read(5,100) achar 100 format(a1) if (achar .eq. "1" ) then print *,' restarting' call rs_read (ddd, time, ncycle, numdmp, 4, nxzones, nyzones, nbdy ) print *,' ncycle=',ncycle print *,' numdmp=',numdmp else c initialize counters ncycle = 0 time = 0.0 numdmp = 0 c c initialize the hydrodynamic quantities: c call set2dshock( ddd, nxzones, nyzones, nbdy ) c endif c c Request the number of timesteps to complete print *,' Please enter number of timesteps to perform.' read *,nstop c c Adjust nstop so that it is a multiple of ndump. This avoids waste. c nstop = ndump * ( (nstop +ndump -1)/ndump ) print *,' Thank you. The program' print *,' will stop after ',nstop,' timesteps.' c c add the number of completed timesteps ncycle nstop = nstop + ncycle tstop = float (nstop) * DT c c CONTROL performs the entire hydrodynamic calculation. c dtime = DT gamma = EOSGAM call cntrol ( xl, yl, ddd, gamma, dtime, smlrho, smalle, smallp, ^ smallu, ncycle, numdmp, time, ndump, tdump, nstop, tstop, ^ nxzones, nyzones, nbdy) c write a restart dump before exiting. call rs_write (ddd, time, ncycle, numdmp, 4, nxzones, nyzones, nbdy ) stop end subroutine rs_write (ddd, time, ncycle, numdmp, nv, nx, ny, nbdy ) c c This subroutine writes a minimal restart dump. A restart dump contains c all the information necessary to restart the program from the point c at which the restart dump was written. c dimension ddd( 1:NV, 1-NBDY:NX+NBDY, 1-NBDY:NY+NBDY ) c c ddd = the 3-D array into which the NV principal 2-D arrays of c basic physical variables are interleaved in order to c optimize memory bandwidth on microprocessor systems. c These are, at a minimum: c ddd(1,i,j) = rrho(i,j) = the zone-averaged densities. c ddd(2,i,j) = pp(i,j) = the zone-averaged pressures. c ddd(3,i,j) = uux(i,j) = the zone-averaged x-velocities. c ddd(4,i,j) = uuy(i,j) = the zone-averaged y-velocities. c c c numdmp = the sequential number of the nextimage 0->whatever, numdmp c NV = number of different zone averaged quantities in DDD 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 ncycle = the number of completed timesteps c time = the program time. c c Open the file "RESTART" iunit = 34 call ppm98_binopen_trunc ("RESTART", iunit) c c write out the scalars first c nbytes = 4 call ppm98_binwrite (iunit, ncycle, nbytes) call ppm98_binwrite (iunit, numdmp, nbytes) call ppm98_binwrite (iunit, time, nbytes) c c now write out ddd c nbytes = 4 * nv * (nx+2*nbdy)* (ny+2*nbdy) call ppm98_binwrite (iunit, ddd, nbytes) call ppm98_binclose (iunit) return end c subroutine rs_read (ddd, time, ncycle, numdmp, nv, nx, ny, nbdy ) c c This subroutine reads a minimal restart dump. A restart dump contains c all the information necessary to restart the program from the point c at which the restart dump was written. c dimension ddd( 1:NV, 1-NBDY:NX+NBDY, 1-NBDY:NY+NBDY ) c c ddd = the 3-D array into which the NV principal 2-D arrays of c basic physical variables are interleaved in order to c optimize memory bandwidth on microprocessor systems. c These are, at a minimum: c ddd(1,i,j) = rrho(i,j) = the zone-averaged densities. c ddd(2,i,j) = pp(i,j) = the zone-averaged pressures. c ddd(3,i,j) = uux(i,j) = the zone-averaged x-velocities. c ddd(4,i,j) = uuy(i,j) = the zone-averaged y-velocities. c c c NV = number of different zone averaged quantities in DDD 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 numdmp = the sequential number of the nextimage 0->whatever, numdmp c ncycle = the number of completed timesteps c time = the program time. c c Open the file "RESTART" iunit = 34 call ppm98_binopen ("RESTART", iunit) c c read the scalars first c nbytes = 4 call ppm98_binread (iunit, ncycle, nbytes) call ppm98_binread (iunit, numdmp, nbytes) call ppm98_binread (iunit, time, nbytes) c c now read ddd c nbytes = 4 * nv * (nx+2*nbdy)* (ny+2*nbdy) call ppm98_binread (iunit, ddd, nbytes) call ppm98_binclose (iunit) return end subroutine cntrol ( xl, yl, ddd, gamma, dtime, smlrho, smalle, ^ smallp, smallu, ncycle, numdmp, time, 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 5 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 numdmp = a count of the number of image dumps completed so far c ncycle = a count of the number of time steps taken so far c time = the elapsed problem time to this point 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 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 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 c Create an image of the initial state. The initial state will be c sequence number 0. if ( numdmp .eq. 0 ) then 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. endif c ntodmp = ndump ttodump = tdump c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c The 1000 loop advances the calculation by two timesteps per c c iteration c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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) return if (time .ge. tstop) return 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 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) = pp(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) 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) = pq1 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) = pq2 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) = pq3 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) = pq4 ddd(3,i,k) = uxq4 ddd(4,i,k) = uyq4 400 continue c c 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), ^ 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,4) 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 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 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) = pp(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 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 character*1 igray(2*nx,2*ny) character*1 acharc, acharc2 character*8 namdmp 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) = pp(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 pp = the zone-averaged pressures. 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 acharc = 'PREFIX' acharc2 = 'e' c c determine individual min/max for scaling c rhomin = ddd(1,1,1) rhomax = ddd(1,1,1) pmin = ddd(2,1,1) pmax = ddd(2,1,1) uxmin = ddd(3,1,1) uxmax = ddd(3,1,1) uymin = ddd(4,1,1) uymax = ddd(4,1,1) do 1000 j = 1,ny do 1000 i = 1,nx ii = ii + 1 rhomin = min( rhomin,ddd(1,i,j)) rhomax = max( rhomax,ddd(1,i,j)) pmin = min( pmin,ddd(2,i,j)) pmax = max( pmax,ddd(2,i,j)) uxmin = min( uxmin,ddd(3,i,j)) uxmax = max( uxmax,ddd(3,i,j)) uymin = min( uymin,ddd(4,i,j)) uymax = max( uymax,ddd(4,i,j)) 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 c irhooff = 0 ipoff = nx iuxoff = 2*nx*ny iuyoff = 2*nx*ny + nx do j = 1,2*ny do i = 1,2*nx igray(i,j) = char(100) enddo enddo ii = 0 do j = 1,ny do i = 1,nx ii = ii + 1 rthyng = rhofac * (ddd(1,i,j) - rhomin) pthyng = pfac * (ddd(2,i,j) - pmin) uxthyng = uxfac * (ddd(3,i,j) - uxmin) uythyng = uyfac * (ddd(3,i,j) - uymin) igray(i,j) = char(int(rthyng)) igray(i+nx,j) = char(int(pthyng)) igray(i,j+ny) = char(int(uxthyng)) igray(i+nx,j+ny) = char(int(uythyng)) enddo enddo c c c c Construct name of image file and store in NAMDMP: c write (namdmp,1004) acharc, acharc2, numdmp 1004 format (a1,a1,'im',i4.4) print *,' namdmp=',namdmp c c c now write it to file namdmp iunit = 34 nbytes = 4*nx*ny call ppm98_binopen_trunc (namdmp, iunit) call ppm98_binwrite (iunit, igray, nbytes) call ppm98_binclose (iunit) c numdmp = numdmp + 1 c c return end