c for single precision: c c m4 -B6666 -Difdouble=0 Driverlr1dgpaqc.EXA > driver.f c c for double precision: c c m4 -B6666 -Difdouble=1 Driverlr1dgpaqc.EXA > driver.f c c======================================================================= c= = c c C WHAT This program is an example of a simple driver to c compute hydrodynamics in lag + remap style c with passively advected quantities (PAQs), using c routines from the PPMLIB library. The inital c problem consists of Sod's Shock tube with reflecting c boundaries. Allowance for example passively advected c quantities is provided. No input is required. Just c compile and run, linking with the PPMLIB library. c c c ASSUMPTIONS 1 dimension c c Cartesian coordinates (planar symmetry) c c Nonuniform 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 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 PAQS Allowance of passively adveced quantities (PAQ) is make. c a PAQ is a quantity which is merely advected by the c flow without direct interaction with the hydrodynamics, c (examples: smoke or dye) or, if interaction exists, c such interaction is separable (example: transveres c velocities.) Both volume weighted and mass weighted c advection is allowed. c c c 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 DIMENSIONS c Three types of arrays exist: zone interface centered, c 1 dimensional strips of zone centered data, and two c dimensional "packed" arrays of 1 dimensional strips of c zone centered data. 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. The final array may be thought of as a c two dimensional array. Most arrays are zone centered. c c It will be assumed that the dimensions of all zone-averaged c type strips are C [1-NBDY:N+NBDY], c and the dimensions of all PAQ arrays are C [1-NBDY:N+NBDY,1:NPAQ], 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 NPAQ is the total number of PAQS in the c calculation. 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 the new PAQS [i.e. PAQNU, will have updated values on the c range C [1:N,1:NPAQ], C and any returned interface type thingies will have valid C values on the range C [1:N+1]. c c c 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), (paq(i,j),j=1,npaq) c (where npaq is the total number of PAQs) 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 program Driverlr1dgpaqc c c This is a simple 1-D triver to solve a Riemann Shock Tube in Lagrangian c + Remap Style. 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) parameter (NBYV = 3) parameter (NBYM = 3) c Note that because of the fixed timestep, N, XMAX, and DT are not c independent. c Note that (NBYV+NBYM) must equal 1 or greater) 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 NBYV = the number of volume weighted PAQs. c NBYM = the number of mass weighted PAQs. c c Do not change this one: parameter (NBDY = 9) c ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) c c dimension xl(1-NBDY:NZONES+NBDY+1) dimension rho(1-NBDY:NZONES+NBDY) dimension p(1-NBDY:NZONES+NBDY) dimension u(1-NBDY:NZONES+NBDY) dimension rhonu(1-NBDY:NZONES+NBDY) dimension pnu(1-NBDY:NZONES+NBDY) dimension unu(1-NBDY:NZONES+NBDY) dimension paq(1-NBDY:NZONES+NBDY,1:(NBYV+NBYM)) dimension paqnu(1-NBDY:NZONES+NBDY,1:(NBYV+NBYM)) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SET constant parameters c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c if ((NBYV+NBYM) .lt. 1) then print *,' The initialization parameters NBYV and NBYM seem' print *,' confused. These parameters indicate the number' print *,' of volume weighted and mass weighted passively' print *,' advected quantities, resp. Their sum should be' print *,' at least 1, otherwise unpredictable results may' print *,' be obtained. Current values are:' print *,' NBYV=',NBYV print *,' NBYM=',NBYM print *,' (NBYV+NBYM)=',NBYV+NBYM print *,' If you wish to continue, enter 1. Any other input' print *,' terminates life.' read *,ilife if( ilife .ne. 1) then call abort() stop endif endif c c For round off ... smlrho = 1.e-06 smalle = 1.e-06 smallp = 1.e-06 smallu = 1.e-06 small = 1.e-06 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SET initial state: c c sod 1-D shock tube c c c p left = 1.0 p right = 0.1 c rho left = 1.0 rho right = 0.125 c u left = 0.0 u right = 0.0 c gamma = 1.4 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc gamma = EOSGAM dtime = DT n = nzones TMAX = float( MAXSTEPS ) * dtime xlimit = XMAX deex = xlimit / float(N) SINAMP = 0.2 SINFAC = 6.28318530717958 / 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 c initialize passively advected quantities. The Ith PAQ is based on c c cos(x) ** I c c nbyvol = NBYV nbymass = NBYM sfac = 6.*3.14159265358979 / float(N) DO IPAQ=1,nbyvol+nbymass DO I=1,N ghoti = cos(float(i)*sfac) paq(i,ipaq) = float(ipaq)*( ghoti ** ipaq ) ENDDO ENDDO c c time = 0.0 nsteps = 0 1010 continue c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Boundary Condition Step c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Impose reflecting boundaries: c DO I=1-NBDY, 0 XL(I) = 2.0 * XL(1) - XL(2 - I) RHO(I) = RHO(1-I) U(I) = -U(1-I) Ctest U(I) = U(1-I) P(I) = P(1-I) do ipaq = 1,nbyvol+nbymass PAQ(I,IPAQ) = PAQ(1-I,IPAQ) enddo ENDDO c 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) Ctest U(I) = U(2*N+1-I) P(I) = P(2*N+1-I) do ipaq = 1,nbyvol+nbymass PAQ(I,IPAQ) = PAQ(2*N+1-I,IPAQ) enddo ENDDO c c courmx = 0.0 c c Compute a single timestep direct eulerian hydro update c nn = N nnbdy = NBDY k = 1-nbdy ibug_bke = 0 call ifelse(ifdouble,1,d_, ) do_ppmlr_1dc_paq_gamma ^ (xl(k), rho(k), p(k), u(k), ^ paq, rhonu(K), pnu(k), unu(k), paqnu, gamma, dtime, ^ smlrho, smallp, smallu, smalle, ^ courmx, nbyvol, nbymass, nn, nnbdy) 666 format (i3, 1h ,1p, 8(1x,e12.5)) c c c Rename for the next timestep print *,' ' print *,' ' print *,' ' print *,' ' print *,' (non-PAQ) Updated results are:' print *,' i xl(i) rho(i) p(i) u(i) ' do 9999 i=1,n rho(i) = rhonu(i) u(i) = unu(i) p(i) = pnu(i) write(*,666) i,xl(i),rho(i),p(i),u(i) 9999 continue do ipaq=1, nbyvol + nbymass do i=1,n paq(i,ipaq) = paqnu(i,ipaq) enddo enddo time = time + dtime nsteps = nsteps+1 print *,' ****** NSTEP TIME COURMX **** = ',nsteps,time,courmx c wreate some output files for mongo , one every timestep, c one every 10 timesteps, one every 25 timesteps: npaq = nbyvol + nbymass ido10 = nsteps - 10*(nsteps/10) ido25 = nsteps - 25*(nsteps/25) 669 format (i3, 1h ,1p, 12(1x,e9.2)) write(101,669) i,xl(1), 0.0, 0.0, 0.0, (0.0,i=1,npaq) if (ido10 .eq. 0) write(110,669) i,xl(1), 0.0, 0.0, 0.0, (0.0,i=1,npaq) if (ido25 .eq. 0) write(125,669) i,xl(1), 0.0, 0.0, 0.0, (0.0,i=1,npaq) do i=1,n write(101,669) i,xl(i),rho(i),p(i),u(i),(paq(i,j),j=1,npaq) if (ido10 .eq. 0) write(110,669) i,xl(i),rho(i),p(i),u(i), (paq(i,j),j=1,npaq) if (ido25 .eq. 0) write(125,669) i,xl(i),rho(i),p(i),u(i), (paq(i,j),j=1,npaq) enddo write(101,669) i,xl(n+1), 0.0, 0.0, 0.0, (0.0,i=1,npaq) if (ido10 .eq. 0) write(110,669) i,xl(n+1), 0.0, 0.0, 0.0, (0.0,i=1,npaq) if (ido25 .eq. 0) write(125,669) i,xl(n+1), 0.0, 0.0, 0.0, (0.0,i=1,npaq) c if(time .lt. tmax) goto 1010 if(nsteps .lt. MAXSTEPS) goto 1010 close(101) close(110) close(125) c || || c || || c || || c || || c ** || || c ** || || c || || c || || c || || c ** || || c ** || || c || || c || || c || || c || || stop end c c c c