c c define(IQX,59) define(IQY,61) define(IQZ,67) define(NBDY,9) define(NVAR,6) define(MAXCOMMENTS,6) program cdumps3dexample c c WHAT: This program is designed to test and illustrate the PPMLIB c compressed dump routines suitable to three-dimensional data sets. c It constructs a data set similar to that found in a typical c PPMLIB hydrodynamic calculation. That is: c a) A skin of fake or boundary zones is appended to all c sides of the data. c b) The computational grid is defined by the locations of the lefthand c zone interfaces, stored in 1-d arrays XINTRL YINTRL, AND ZINTRL. c In order to fully define the computational grid using lefthand zone c interfaces, one addional interface is required on the right of c these arrays. They are therefore dimensioned one greater on the c right than would otherwise be expected. c c) The three-dimensional arrays are stored in interleaved order in a c four-dimensional array DDD. This interleaved ordering provides c for more efficient memory access on cache based machines. c It is assumed that the ordering of the values in DDD is: density c pressure, x-velocity, y-velocity, z-velocity, diagnostic. Thus c the array elements for zone DDD(1,i,j), DDD(2,i,j), DDD(3,i,j), c DDD(4,i,j), DDD(5,i,j), and DDD(6,i,j) refer to the values of c density, pressure, x-velocity, y-velocity, z-velocity, and the c diagnostic for zone (i,j), respectively. This specific ordering c is arbitrary and determined by the programmer. c c SPECIAL COMMENTS: It is assumed that this file will be preprocessed using c the m4 macro processor with the option -Difdouble=0 for single c precision or -Difdouble=1 for double precision. That is: c single precision: m4 -Difdouble=0 cdumps3dexample.m4 > cdumps3dexample.f c double precision: m4 -Difdouble=1 cdumps3dexample.m4 > cdumps3dexample.f ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) c ifelse(ifdouble,1,8,4) parameter(nxzones=IQX) parameter(nyzones=IQY) parameter(nzzones=IQZ) parameter(nbdy=NBDY) parameter(nvar=NVAR) parameter(maxcomments=MAXCOMMENTS) c c These are 1d arrays which contain information defining the grid, c that is, the locations of the left hand zone interfaces. One c extra zone is requires on the left to fully specify the grid. dimension xintrl(1-nbdy:nxzones+nbdy+1) dimension yintrl(1-nbdy:nyzones+nbdy+1) dimension zintrl(1-nbdy:nzzones+nbdy+1) c character*20 filename c c This will hold some fascinating commentary describing the problem c for inclusion in the compressed dump file. character*80 comments(maxcomments) c c This array contains zone averaged information of the state variables c in interlieaved order, making it a 4d array for three spatial dimensions c They are: c ddd(1,*,*,*) = density c ddd(2,*,*,*) = pressure c ddd(3,*,*,*) = x-velocity c ddd(4,*,*,*) = y-velocity c ddd(5,*,*,*) = z-velocity c ddd(6,*,*,*) = diagnostic (an arbitrary quantity) dimension ddd(nvar,1-nbdy:nxzones+nbdy, ^ 1-nbdy:nyzones+nbdy, ^ 1-nbdy:nzzones+nbdy) c c Initialize the disk file name filename = 'p14a0000' c c Initialize the commentary comments(1)=' This is the result from the program cdumps3dexample' comments(2)=' Its purpose is to demonstrate the PPMLIB compressed' comments(3)=' dumps routines. The arrays are REAL*ifelse(ifdouble,1,8,4)' ncomments = 3 c c Initialize the grid call init_grid(xintrl, yintrl, zintrl, nxzones, nyzones, nzzones, nbdy) c c Initialize the zone average grid call init_ddd(ddd,nvar,nxzones,nyzones,nzzones,nbdy) c c Open the file, write the comments including the date of file c creation. iunit = 66 call ppm98_binopen_trunc(filename,iunit) do i=1,ncomments call ppm98_write_comment(comments(i),iunit) enddo call ppm98_timestamp(iunit) c This illustrates adding a constant value to the compressed dump file gamma= 1.4ifelse(ifdouble,1,d,e)+00 call ifelse(ifdouble,1,d_,)ppm98_var0d('ratio of specific heats', ^ 'GAMMA',gamma,iunit) c c Now record the grid information in the compressed dump routine. The c three axes will be assigned the names 'X', 'Y', and 'Z'. These c should be capitol letters, as their derivatives will be references c with the equivalent lower case letters in a3d. call ifelse(ifdouble,1,d_,)ppm98_mesho3('X', 'Y', 'Z', nbdy, nxzones, nbdy, ^ nbdy, nyzones, nbdy, ^ nbdy, nzzones, nbdy, ^ xintrl, yintrl, zintrl, iunit ) c c Now add the individual zone average quantities to the compressed dump call ifelse(ifdouble,1,d_,)ppm98_var3d('Density', 'RHO', ^ nbdy, nxzones, nbdy, ^ nbdy, nyzones, nbdy, ^ nbdy, nzzones, nbdy, ^ nvar, 1, ddd, 0, iunit) call ifelse(ifdouble,1,d_,)ppm98_var3d('Pressure', 'P', ^ nbdy, nxzones, nbdy, ^ nbdy, nyzones, nbdy, ^ nbdy, nzzones, nbdy, ^ nvar, 2, ddd, 0, iunit) call ifelse(ifdouble,1,d_,)ppm98_var3d('X velocity', 'VX', ^ nbdy, nxzones, nbdy, ^ nbdy, nyzones, nbdy, ^ nbdy, nzzones, nbdy, ^ nvar, 3, ddd, 0, iunit) call ifelse(ifdouble,1,d_,)ppm98_var3d('Y velocity', 'VY', ^ nbdy, nxzones, nbdy, ^ nbdy, nyzones, nbdy, ^ nbdy, nzzones, nbdy, ^ nvar, 4, ddd, 0, iunit) call ifelse(ifdouble,1,d_,)ppm98_var3d('Z velocity', 'VZ', ^ nbdy, nxzones, nbdy, ^ nbdy, nyzones, nbdy, ^ nbdy, nzzones, nbdy, ^ nvar, 5, ddd, 0, iunit) call ifelse(ifdouble,1,d_,)ppm98_var3d( ^ '3D chess board diagnostic', 'DIAG', ^ nbdy, nxzones, nbdy, ^ nbdy, nyzones, nbdy, ^ nbdy, nzzones, nbdy, ^ nvar, 6, ddd, 0, iunit) c Close the file and exit call ppm98_binclose(iunit) end subroutine init_ddd(ddd,nvar,nx,ny,nz,nbdy) ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension ddd(nvar, 1-nbdy:nx+nbdy, ^ 1-nbdy:ny+nbdy, ^ 1-nbdy:nz+nbdy) c c ddd(1,*,*,*) = density c ddd(2,*,*,*) = pressure c ddd(3,*,*,*) = x-velocity c ddd(4,*,*,*) = y-velocity c ddd(5,*,*,*) = z-velocity c ddd(6,*,*,*) = a diagnostic c c The domain will be divided into octants for density and pressure c c The velocities shall consist of a sine wave twopi = 2.0*atan2(0.,-1.) xfac = twopi/float(nx) yfac = twopi/float(ny) zfac = twopi/float(nz) c c The diagnostic shall consist of alternating cubes of 0.,1. in c a 3d chess boeard pattern. lcube = 16 c xcube = 0.0 do iz=1-nbdy,nz+nbdy if (iz.le.nz/2) then z = 1.0 else z = 2.0 endif zvel = sin(float(iz)*zfac) if (mod(iz,lcube) .eq. 0) xcube = 1.0 - xcube ycube = xcube do iy=1-nbdy,ny+nbdy if (iy.le.ny/2) then y = 3.0 else y = 4.0 endif yvel = sin(float(iy)*yfac) if (mod(iy,lcube) .eq. 0) ycube = 1.0 - ycube zcube = ycube do ix=1-nbdy,nx+nbdy if (ix.le.nx/2) then x = 5.0 else x = 6.0 endif xvel = sin(float(ix)*xfac) if (mod(ix,lcube) .eq. 0) zcube = 1.0 - zcube ddd(1,ix,iy,iz) = x*y*z ddd(2,ix,iy,iz) = x*y*z ddd(3,ix,iy,iz) = xvel ddd(4,ix,iy,iz) = yvel ddd(5,ix,iy,iz) = zvel ddd(6,ix,iy,iz) = zcube enddo enddo enddo return end subroutine init_grid(x1intrl,x2intrl,x3intrl,n1zones,n2zones,n3zones,nbdy) ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension x1intrl(1-nbdy:n1zones+nbdy+1) dimension x2intrl(1-nbdy:n2zones+nbdy+1) dimension x3intrl(1-nbdy:n3zones+nbdy+1) c Set up a uniform 3d grid where the non boundary portion c ranges from 0 to 1. dx = 1.0 / float(n1zones) do i=1-nbdy,n1zones+nbdy+1 x1intrl(i) = float(i-1)*dx enddo dx = 1.0 / float(n2zones) do i=1-nbdy,n2zones+nbdy+1 x2intrl(i) = float(i-1)*dx enddo dx = 1.0 / float(n3zones) do i=1-nbdy,n3zones+nbdy+1 x3intrl(i) = float(i-1)*dx enddo return end