c c This file is meant ot be preprocessed by "m4" c c for single precision: c c m4 -B6666 -Difdouble=0 Driverdirecte1dgc.EXA > driver.f c c for double precision: c c m4 -B6666 -Difdouble=1 Driverdirecte1dgc.EXA > driver.f c======================================================================= c= = c c ************************************************************************** c * THIS CODE IS Driverdirecte1dgc * 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. The inital c problem consists of Sod's Shock tube with reflecting c boundaries. No input is required. Just compile and c run, linking with the PPMLIB library. c c ASSUMPTIONS 1 dimension c c Cartesian coordinates (planar symmetry) 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 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 DIMENSIONS c Two types of 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 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], C and 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 of an ASCII text in which each c line consists of c i, xl(i), rho(i), p(I), u(I) c displayed on the user's screen after each timestep c and in the files fort.101, fort.110, fort.125, for c every timestep, every 10 timesteps, and every 25 c timesteps. The file fort.999 contains the results c from the last timestep only. The format of the fort.* c is compatable with the MONGO plotting package popular c in the astronomical community. c C program Driverdirecte1dg0c c c c some ability to change the program parameters is provided via: c parameter (NZONES = 128) parameter (MAXSTEPS = 100) parameter (EOSGAM = 1.4) parameter (DT = 0.0025) parameter (XMAX = 1.0) c Note that because of the fixed timestep, N, XMAX, and DT are not c independent. c c NZONES = number of real computational zones in the calculation 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 value of the grid. The grid runs from (0,XMAX)c c c Do not change this one: parameter (NBDY = 5) c c dimension xl(1-NBDY:NZONES+NBDY+1) dimension rho(1-NBDY:NZONES+NBDY), p(1-NBDY:NZONES+NBDY) dimension u(1-NBDY:NZONES+NBDY) dimension rhonu(1-NBDY:NZONES+NBDY), unu(1-NBDY:NZONES+NBDY) dimension pnu(1-NBDY:NZONES+NBDY) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SET constant parameters c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c c For round off ... smlrho = 1.e-06 smalle = 1.e-06 smallp = 1.e-06 smallu = 1.e-06 gamma = EOSGAM dtime = DT xlimit = XMAX TMAX = float( MAXSTEPS ) * dtime n = nzones cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c SETUP initial state: c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc k = 1-nbdy call setup( xl(k), rho(k), p(k), u(k), xlimit, gamma, n, nbdy ) c c Initialize counters for total time and number of timestep: time = 0.0 nsteps = 0 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c The 1010 loop is the hydro loop. One iteration is one timestep c c Each timestep is divided into three substeps: c c 1: Impose boundary conditions c c 2: Call a PPN hydro package to produce new zone averages c c 3: Replace old zone averages with new in preparation c c for the next timestep. c c 4: Produce output if desired c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 1010 continue c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Step 1: Boundary Condition Step c c Boundary conditions are set at the beginning of each timestep c c by placing appropriate values of XL, RHO, P, and U in the fake c c zones on the left and right ends of the arrays. c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Impose reflecting boundaries: k=1-nbdy call boundary_refl(xl(k), rho(k), u(k), p(k), n, nbdy) c c clear old maximum Courant number: courmx = 0.0 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Step 2: Perform a single timestep direct eulerian hydro update c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c nn = N nnbdy = NBDY k = 1-nbdy call do_ppmde_1dc_gamma (xl(k), rho(k), p(k), u(k), rhonu(k), ^ pnu(k), unu(k), gamma, dtime, smlrho, smallp, smallu, smalle, ^ courmx, nn, nnbdy) c c Increment counters: time = time + dtime nsteps = nsteps+1 c Rename for the next timestep (It is assumed that XL is unchanged) c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Step 3: Replace old zone averages with new zone averages c c i preparation of the next timestep. c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do i=1,n rho(i) = rhonu(i) u(i) = unu(i) p(i) = pnu(i) enddo cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Step 4: Output c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c A: print values of updated quantities to screed each timestep print *,' ' print *,' ' print *,' ' print *,' Updated results are:' print *,' i xl(i) rho(i) p(i) u(i) ' do 9999 i=1,n write(*,666) i,xl(i),rho(i),p(i),u(i) 666 format (i3, 1h ,1p, 8(1x,e12.5)) 9999 continue print *,' ****** NSTEP TIME COURMX **** = ',nsteps,time,courmx c B: Create disk data files: c create some output files for mongo , one every timestep, c one every 10 timesteps, one every 25 timesteps: ido10 = nsteps - 10*(nsteps/10) ido25 = nsteps - 25*(nsteps/25) write(101,666) i,xl(1), 0.0, 0.0, 0.0 if (ido10 .eq. 0) write(110,666) i,xl(1), 0.0, 0.0, 0.0 if (ido25 .eq. 0) write(125,666) i,xl(1), 0.0, 0.0, 0.0 do i=1,n write(101,666) i,xl(i),rho(i),p(i),u(i) if (ido10 .eq. 0) write(110,666) i,xl(i),rho(i),p(i),u(i) if (ido25 .eq. 0) write(125,666) i,xl(i),rho(i),p(i),u(i) enddo write(101,666) i,xl(n+1), 0.0, 0.0, 0.0 if (ido10 .eq. 0) write(110,666) i,xl(n+1), 0.0, 0.0, 0.0 if (ido25 .eq. 0) write(125,666) i,xl(n+1), 0.0, 0.0, 0.0 c Go back to 1010 if further timesteps are to be computed: c c if(time .lt. tmax) goto 1010 if(nsteps .lt. MAXSTEPS) goto 1010 c close the disk files and stop close(101) close(110) close(125) c stop end c c c c cSubroutineSubroutineSubroutineSubroutineSubroutineSubroutineSubroutineSubroutinc c #### # # ##### ##### #### # # ##### # # # ###### c # # # # # # # # # # # # # ## # # c #### # # ##### # # # # # # # # # # # ##### c # # # # # ##### # # # # # # # # # # c # # # # # # # # # # # # # # # ## # c #### #### ##### # # #### #### # # # # ###### c cSubroutineSubroutineSubroutineSubroutineSubroutineSubroutineSubroutineSubroutinc c subroutine setup( xl, rho, p, u, xmax, gamma, n, nbdy ) c c This routine sets up a 1d shock tube problem on a nonuniform grid. Said c grid ranges from xl=0 to xl=xmax. The zone widths have a smooth sinusoidal c variation. c c xl = locations of lefthand zone interfaces c rho = zone averaged density c p = zone averaged pressure c u = zone averaged velocity c c n = number of real zones in calculation c nbdy = number of boundary zones added to each end of 1d arrays to c implement boundary conditions. c dimension xl(1-nbdy:n+nbdy+1) dimension rho(1-nbdy:n+nbdy) dimension p(1-nbdy:n+nbdy) dimension u(1-nbdy:n+nbdy) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c SET initial state: c c sod 1-D shock tube c c p left = 1.0 p right = 0.1 c c rho left = 1.0 rho right = 0.125 c c u left = 0.0 u right = 0.0 c c gammaa = 1.4 c c c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc sinamp = 0.2 sinfac = 6.28318530717958 / float (N) deex = xmax / float(N) xl(1) = 0.0 DO I=2,N+1 xl(i) = xl(i-1) + deex * ( 1 + sinamp * sin(sinfac*float(i)) ) ENDDO DO I=1,N IF ( I .LE. N/2 ) THEN P(I) = 1.0 RHO(I) = 1.0 U(I) = 0.0 ELSE P(I) = 0.1 RHO(I) = 0.125 U(I) = 0.0 ENDIF c pnu(i) = p(i) c print *,'irpue=',i,rho(i),p(I),u(I) ENDDO c return end c c c c cSubroutineSubroutineSubroutineSubroutineSubroutineSubroutineSubroutineSubroutinc c #### # # ##### ##### #### # # ##### # # # ###### c # # # # # # # # # # # # # ## # # c #### # # ##### # # # # # # # # # # # ##### c # # # # # ##### # # # # # # # # # # c # # # # # # # # # # # # # # # ## # c #### #### ##### # # #### #### # # # # ###### c cSubroutineSubroutineSubroutineSubroutineSubroutineSubroutineSubroutineSubroutinc subroutine boundary_refl(xl, rho, u, p, n, nbdy) c c This routine imposed reflecting boundary conditions on c 1d strips c c xl = locations of lefthand zone interfaces c rho = zone averaged density c p = zone averaged pressure c u = zone averaged velocity c c n = number of real zones in calculation c nbdy = number of boundary zones added to each end of 1d arrays to c implement boundary conditions. c dimension xl(1-nbdy:n+nbdy+1) dimension rho(1-nbdy:n+nbdy) dimension p(1-nbdy:n+nbdy) dimension u(1-nbdy:n+nbdy) c Impose reflecting boundaries: c c Left: DO I=1-NBDY, 0 XL(I) = 2.0 * XL(1) - XL(2 - I) RHO(I) = RHO(1-I) U(I) = -U(1-I) P(I) = P(1-I) ENDDO c c c Right: DO I = N+1, N+NBDY XL(I+1) = 2.0 * XL(N+1) - XL(2*N+1 - I) RHO(I) = RHO(2*N+1-I) U(I) = -U(2*N+1-I) P(I) = P(2*N+1-I) ENDDO return end c