program ppm implicit none c..one-dimensional explicit hydrodynamic code using the c..lagrange + remap formulation of the piecewise-parabolic method. c..programmed by bruce fryxell, december 1992 c..cleaned by frank timmes march 2005 c..bring in the common blocks include 'ppm.dek' c..set up grid and initial distribution of conserved quantities call grid call init call tstep write(6,*) nr c..main time loop do nstep = 1, nmax write (6,*) nstep, t, dt if (t + dt .gt. tmax) dt = tmax - t t = t + dt call get call hydro call store call tstep if (t .ge. tmax) go to 10 enddo 10 continue c..write out results call output stop end subroutine grid implicit none save c..compute grid coordinates c..bring in the common blocks include 'ppm.dek' c..local variables integer i double precision grdlnr,delr,delr2 c..choose coordinate system c..alpha = 0 ==> cartesian c..alphar = 1 ==> cylindrical alpha = 0 c..choose the number of grid points nr = 500 c..set up grid in r (or x) direction c..the following code assumes a uniformly spaced grid. c..in practice, any grid spacing may be used. c..here we set up a grid between zero and one. grdlnr = 5.0d0 delr = grdlnr / nr delr2 = 0.5d0 * delr do i = 1, nr rr(i) = i * delr rl(i) = rr(i) - delr r (i) = rl(i) + delr2 c write(6,115) i,rl(i),r(i),rr(i) c 115 format(1x,i5,1p3e12.4) enddo c read(5,*) return end subroutine init implicit none save c..define distributions of conserved quantities at the initial time. c..bring in the common blocks include 'ppm.dek' c..local variables integer i,njmpl,njmpr double precision rorite,roleft,urite,uleft,pleft,prite double precision pi parameter (pi = 3.1415926535898d0) c..the coding below sets up the initial distribution for c..calculating a rayleigh-taylor instability between two planar c..layers of fluid in a constant gravitational field. c..define some initial constants c.. nmax = number of time steps allowed before code stops c.. t = initial physical time c.. dt = guess for initial time step value c.. dtmin, dtmax = minimum and maximum values allowed for time step c.. cfl = fraction of maximum stable time step to use c.. nriem = maximum number of iterations in riemann solution c.. gamma = the ratio of specific heats nmax = 1000000 t = 0.0d0 tmax = 0.40d0 dt = 1.0d-3 dtmin = 1.d-10 dtmax = tmax cfl = 0.8d0 nriem = 10 gamma = 1.4d0 c..the following constants are 'small' values for dimensionless c..numbers, density, pressure, and energy. they should c..be set to values a few orders of magnitude less than values c..expected for the corresponding variable. they are used primarily c..to prevent errors such as dividing by zero which can result c..if a variable undershoots its correct value. small = 1.d-10 smlrho = 1.d-10 smallp = 1.d-10 smalle = 1.d-10 gamm1 = gamma - 1.d0 gamp1 = gamma + 1.d0 gmfc = 0.5d0 * gamm1 gamfac = 0.5d0 * gamp1 / gamma c..gr is the gravitational accelerations gr = 0.d0 c..define the values of the conserved quantities at each grid point pleft = 100.0d0 uleft = 0.0d0 roleft = 10.d0 prite = 1.0d0 urite = 0.0d0 rorite = 1.0d0 njmpl = nr / 2 njmpr = nr - njmpl do i = 1, nr if ( r(i) .le. 2.0d0) then densty(i) = roleft vr(i) = uleft press(i) = pleft energy(i) = press(i) / (gamm1 * densty(i)) else densty(i) = rorite vr(i) = urite press(i) = prite energy(i) = press(i) / (gamm1 * densty(i)) endif enddo return end subroutine get implicit none save c..copy a row of the grid from the three-dimensional arrays c..into the one-dimensional arrays c..bring in the common blocks include 'ppm.dek' c..local variables integer i,i4,i5,i9 c..initialize n = nr np1 = n + 1 np2 = n + 2 np3 = n + 3 np4 = n + 4 np5 = n + 5 np6 = n + 6 np7 = n + 7 np8 = n + 8 g = gr do i = 1, nr i4 = i + 4 rho(i4) = densty(i) u(i4) = vr(i) e(i4) = energy(i) p(i4) = press(i) xl(i4) = rl(i) x(i4) = r(i) xr(i4) = rr(i) dx(i4) = xr(i4) - xl(i4) enddo c..reflecting boundary conditions do i = 1, 4 i9 = 9 - i rho(i) = rho(i9) u(i) = -u(i9) e(i) = e(i9) p(i) = p(i9) dx(i) = dx(i9) enddo do i = 1, 4 i5 = 5 - i xl(i5) = xl(i5+1) - dx(i5) xr(i5) = xr(i5+1) - dx(i5+1) x(i5) = 0.5d0 * (xl(i5) + xr(i5)) enddo do i = np5, np8 i9 = 2 * n + 9 - i rho(i) = rho(i9) u(i) = -u(i9) e(i) = e(i9) p(i) = p(i9) dx(i) = dx(i9) enddo do i = np5, np8 xl(i) = xl(i-1) + dx(i-1) xr(i) = xr(i-1) + dx(i) x(i) = 0.5d0 * (xl(i) + xr(i)) enddo c.. do i = 1, np8 tau(i) = 1.d00 / rho(i) c(i) = sqrt (gamma * p(i) * rho(i)) dm(i) = rho(i) * dx(i) enddo if (alpha .eq. 1) then do i = 1, np8 dm(i) = dm(i) * dabs(x(i)) enddo endif return end subroutine store implicit none save c..store the updated values of the one-dimensional arrays into c..the appropriate row of the three-dimensional arrays c..bring in the common blocks include 'ppm.dek' c..local variables integer i,i4 do i = 1, nr i4 = i + 4 densty(i) = rhonu(i4) vr(i) = unu(i4) energy(i) = enu(i4) press(i) = pnu(i4) rl(i) = xlnu(i4) r(i) = xnu(i4) rr(i) = xrnu(i4) enddo return end subroutine hydro implicit none save c..solve the lagrangian hydrodynamic equations on c..a single row or column of the grid c..bring in the common blocks include 'ppm.dek' c..local variables integer i,imap double precision s1(ndp8), dvol(ndp8), al(ndp8), ar(ndp8), ei do i = 1, np8 dtdx(i) = dt / dx(i) dtdm(i) = dt / dm(i) ce(i) = c(i) / rho(i) enddo imap = 0 call intrfc (imap) call states call rieman do i = 5, np5 xlnu(i) = xl(i) + dt * ustar(i) xrnu(i-1) = xlnu(i) upstar(i) = ustar(i) * pstar(i) enddo do i = 5, np4 dxnu(i) = xrnu(i) - xlnu(i) xnu(i) = 0.5 * (xrnu(i) + xlnu(i)) enddo c..geometry factors if (alpha .eq. 0) then do i = 5, np4 dvol(i) = dxnu(i) al(i) = 1.d0 ar(i) = 1.d0 enddo else do i = 5, np4 dvol(i) = 0.5d0 * dabs (xlnu(i) + xrnu(i)) * dxnu(i) al(i) = 0.5d0 * (xlnu(i) + xl(i)) ar(i) = 0.5d0 * (xrnu(i) + xr(i)) enddo endif do i = 5, np4 rhonu(i) = dm(i) / dvol(i) taunu(i) = 1.d0 / rhonu(i) unu(i) = u(i) - dtdm(i) * 0.5d0 * (al(i) + ar(i)) 1 * (pstar(i+1) - pstar(i)) + dt * g enu(i) = e(i) - dtdm(i) 1 * (ar(i) * upstar(i+1) - al(i) * upstar(i)) 2 + 0.5d0 * dt * (u(i) + unu(i)) * g ei = enu(i) - 0.5d0 * unu(i)**2 ei = max (ei, smalle) pnu(i) = gamm1 * rhonu(i) * ei enddo return end subroutine intrfc (imap) implicit none save c..calculate interface values of density, pressure, and velocity c..bring in the common blocks include 'ppm.dek' c..local variables integer imap,i c..calculate interpolation coefficients call coeff c..perform interpolation call interp (rhol, rho, rhor) call interp ( ul, u, ur) call interp ( pl, p, pr) c..in remap step only, interpolate transverse velocity and c..steepen contact discontinuities if (imap .eq. 1) then call detect endif c..apply monotonicity constraints on interpolation parabolae call monot (rhol, rho, rhor, drho, rho6) call monot ( ul, u, ur, du, u6) call monot ( pl, p, pr, dp, p6) c..flatten interpolation parabolae near shocks which are too thin call flaten c..load the solution arrays do i = 4, np5 rhol(i) = fshk(i) * rho(i) + fshk1(i) * rhol(i) rhor(i) = fshk(i) * rho(i) + fshk1(i) * rhor(i) ul(i) = fshk(i) * u(i) + fshk1(i) * ul(i) ur(i) = fshk(i) * u(i) + fshk1(i) * ur(i) pl(i) = fshk(i) * p(i) + fshk1(i) * pl(i) pr(i) = fshk(i) * p(i) + fshk1(i) * pr(i) enddo do i = 4, np5 drho(i) = rhor(i) - rhol(i) du(i) = ur(i) - ul(i) dp(i) = pr(i) - pl(i) rho6(i) = 6.d0 * rho(i) - 3.d0 * (rhol(i) + rhor(i)) u6(i) = 6.d0 * u(i) - 3.d0 * (ul(i) + ur(i)) p6(i) = 6.d0 * p(i) - 3.d0 * (pl(i) + pr(i)) enddo return end subroutine coeff implicit none save c..calculate coefficients of interpolation polynomial c..bring in the common blocks include 'ppm.dek' c..local variables integer i double precision s1(ndp8), s2(ndp8), s3(ndp8), s4(ndp8) do i = 2, np8 s1(i) = dx(i) + dx(i-1) s2(i) = s1(i) + dx(i) s3(i) = s1(i) + dx(i-1) enddo do i = 2, np7 s4(i) = dx(i) / (s1(i) + dx(i+1)) c1(i) = s4(i) * s3(i) / s1(i+1) c2(i) = s4(i) * s2(i+1) / s1(i) enddo do i = 2, np6 s4(i) = 1.d0 / (s1(i) + s1(i+2)) c3(i) = -s4(i) * dx(i) * s1(i) / s3(i+1) c4(i) = s4(i) * dx(i+1) * s1(i+2) / s2(i+1) c5(i) = (dx(i) - 2.d0 * (dx(i+1) * c3(i) + dx(i) * c4(i))) 1 / s1(i+1) enddo return end subroutine interp (al, a, ar) implicit none save c..interpolate interface values and monotonize c..bring in the common blocks include 'ppm.dek' c..declare the pass double precision al(ndp8), a(ndp8), ar(ndp8) c..local variables integer i double precision s1(ndp8), s2(ndp8), s3(ndp8), armax, armin do i = 2, np8 s1(i) = a(i) - a(i-1) s2(i) = 2.d0 * dabs (s1(i)) enddo do i = 2, np7 dela(i) = c1(i) * s1(i+1) + c2(i) * s1(i) dela(i) = min (dabs (dela(i)), s2(i), s2(i+1)) 1 * sign (1.d0, dela(i)) if (s1(i) * s1(i+1) .le. 0.d0) dela(i) = 0.d0 enddo do i = 2, np6 ar(i) = a(i) + c5(i) * s1(i+1) + c3(i) * dela(i+1) 1 + c4(i) * dela(i) enddo do i = 2, np6 armin = min (a(i), a(i+1)) armax = max (a(i), a(i+1)) ar(i) = min (ar(i), armax) ar(i) = max (ar(i), armin) al(i+1) = ar(i) enddo return end subroutine monot (al, a, ar, da, a6) implicit none save c..apply monotonicity constraint to interpolation parabolae c..bring in the common blocks include 'ppm.dek' c..declare the pass double precision al(ndp8), a(ndp8), ar(ndp8), da(ndp8), a6(ndp8) c..local variables integer i double precision s1,s2,a3 do i = 4, np5 if ((ar(i) - a(i)) * (a(i) - al(i)) .le. 0.d0) then al(i) = a(i) ar(i) = a(i) endif da(i) = ar(i) - al(i) a3 = 3.d0 * a(i) s1 = a3 - 2.d0 * ar(i) s2 = a3 - 2.d0 * al(i) if (da(i) * (al(i) - s1 ) .lt. 0.d0) al(i) = s1 if (da(i) * (s2 - ar(i)) .lt. 0.d0) ar(i) = s2 da(i) = ar(i) - al(i) a6(i) = 6.d0 * a(i) - 3.d0 * (al(i) + ar(i)) enddo return end subroutine flaten implicit none save c..flaten zone structure in regions where shocks are too thin c..bring in the common blocks include 'ppm.dek' c..local variables integer i double precision divu(ndp8), f(ndp8), delp1(ndp8), 1 dpf,dmch0,dmch,s0,sp,delu1,delp2 parameter (dpf = 0.3333333333333333d0) dmch0 = 2.d0 * dpf**2 / (gamma * (2.d0 * gamma + gamp1 * dpf)) dmch0 = sqrt (dmch0) c.. note: the expression for divu in the next loop is valid only for c.. Cartesian coordinates. this should be replaced by the appropriate c.. form for cylindrical coordinates if alpha = 1. however, for c.. most problems, the Cartesian form will provide acceptable results. do i = 3, np6 delu1 = u(i+1) - u(i-1) delp1(i) = p(i+1) - p(i-1) delp2 = p(i+2) - p(i-2) divu(i) = delu1 / (x(i+1) - x(i-1)) sp = abs (delp1(i) / max (delp2, smallp)) s0 = max (0.d0, min (0.5d0, 5.d0 * (sp - 0.75d0))) dmch = abs (delu1 / min (ce(i+1), ce(i-1))) if (dmch .gt. dmch0 .and. divu(i) .lt. 0.d0) then f(i) = s0 else f(i) = 0.d0 endif enddo do i = 4, np5 if (delp1(i) .lt. 0.d0) then fshk(i) = max (f(i), f(i-1)) else fshk(i) = max (f(i), f(i+1)) endif fshk1(i) = 1.d0 - fshk(i) enddo return end subroutine detect implicit none save c..search for contact discontinuities and c..steepen the zone structure if necessary c..bring in the common blocks include 'ppm.dek' c..local variables integer i double precision s1(ndp8), s2(ndp8), s3(ndp8), s4(ndp8), 1 d2rho(ndp8), eta(ndp8), rhodr,rhodl, 2 eta1, eta2, epsln, k0 parameter (eta1 = 20.0d0, eta2 = 0.05d0, 1 epsln = 0.01d0, k0 = 0.1d0) c..load do i = 2, np7 s1(i) = dx(i) + dx(i-1) s2(i) = s1(i) + dx(i+1) s3(i) = rho(i) - rho(i-1) s1(i) = s3(i) / s1(i) enddo do i = 2, np6 d2rho(i) = (s1(i+1) - s1(i)) / s2(i) enddo do i = 2, np8 s1(i) = x(i) - x(i-1) s1(i) = s1(i) * s1(i) * s1(i) enddo do i = 4, np5 s3(i) = - (d2rho(i+1) - d2rho(i-1)) * (s1(i) + s1(i+1)) if (rho(i+1) - rho(i-1) .eq. 0.d0) then s4(i) = small * smlrho else s4(i) = rho(i+1) - rho(i-1) endif eta(i) = s3(i) / ((x(i+1) - x(i-1)) * s4(i)) if (d2rho(i-1) * d2rho(i+1) .gt. 0.d0) eta(i) = 0.d0 enddo do i = 4, np5 s4(i) = epsln * min (dabs (rho(i+1)), dabs (rho(i-1))) 1 - dabs (rho(i+1) - rho(i-1)) if (s4(i) .ge. 0.d0) eta(i) = 0.d0 eta(i) = max (0.d0, min (eta1 * (eta(i) - eta2), 1.d0)) s1(i) = abs (p(i+1) - p(i-1) ) 1 / min (p(i+1), p(i-1) ) s2(i) = abs (rho(i+1) - rho(i-1)) 1 / min (rho(i+1), rho(i-1)) if (gamma * k0 * s2(i) .lt. s1(i)) eta(i) = 0.d0 enddo do i = 4, np5 rhodl = rho(i-1) + 0.5 * dela(i-1) rhodr = rho(i+1) - 0.5 * dela(i+1) rhol(i) = rhol(i) * (1. - eta(i)) + rhodl * eta(i) rhor(i) = rhor(i) * (1. - eta(i)) + rhodr * eta(i) enddo return end subroutine states implicit none save c..compute left and right states for input to riemann problem c..bring in the common blocks include 'ppm.dek' c..local variables integer i double precision s1(ndp8), s2(ndp8), apls, amns, a, forthd parameter (forthd = 1.33333333333333d0) c..load do i = 1, np8 s1(i) = 0.5d0 * min (1.d0, ce(i) * dtdx(i)) s2(i) = 1.d0 - forthd * s1(i) enddo c..geometry factors if (alpha .eq. 1) then do i = 5, np5 apls = xl(i) - s1(i-1) * dx(i-1) a = xl(i) upls(i) = ur(i-1) + (apls * upls(i) - a * ur(i-1)) 1 / (0.5d0 * (apls + a)) enddo else do i = 5, np5 upls(i) = ur(i-1) - s1(i-1) * (du(i-1) - s2(i-1) * u6(i-1)) enddo endif c.. do i = 5, np5 upls(i) = upls(i) + 0.5d00 * dt * g taupls(i) = rhor(i-1) 1 - s1(i-1) * (drho(i-1) - s2(i-1) * rho6(i-1)) taupls(i) = 1.d0 / taupls(i) ppls(i) = pr(i-1) - s1(i-1) * (dp(i-1) - s2(i-1) * p6(i-1)) ppls(i) = max (smallp, ppls(i)) cpls(i) = sqrt (gamma * ppls(i) / taupls(i)) enddo c..geometry factors if (alpha .eq. 1) then do i = 5, np5 amns = xl(i) + s1(i) * dx(i) a = xl(i) umns(i) = ul(i) + (amns * umns(i) - a * ul(i)) 1 / (0.5d0 * (amns + a)) enddo else do i = 5, np5 umns(i) = ul(i) + s1(i) * (du(i) + s2(i) * u6(i)) enddo endif do i = 5, np5 umns(i) = umns(i) + 0.5d00 * dt * g taumns(i) = rhol(i) + s1(i) * (drho(i) + s2(i) * rho6(i)) taumns(i) = 1.d0 / taumns(i) pmns(i) = pl(i) + s1(i) * (dp(i) + s2(i) * p6(i)) pmns(i) = max (smallp, pmns(i)) cmns(i) = sqrt (gamma * pmns(i) / taumns(i)) enddo return end subroutine rieman implicit none save c..solve riemann shock tube problem c..bring in the common blocks include 'ppm.dek' c..local variables integer i,imax, nnn double precision pstr1(ndp8), ustarm, ustarp, emax, zp, zm, 1 error, s1, s2, tol parameter (tol = 1.0d-4) c..first guess (wpls = cpls , wmns = cmns) do i = 5, np5 pstar(i) = pmns(i) - ppls(i) - cmns(i) * (umns(i) - upls(i)) pstar(i) = ppls(i) + pstar(i) * cpls(i) / (cpls(i) + cmns(i)) pstar(i) = max (smallp, pstar(i)) enddo c..begin iteration loop do nnn = 1, nriem do i = 5, np5 s1 = 1.d0 + gamfac * (pstar(i) - ppls(i)) / ppls(i) s2 = 1.d0 + gamfac * (pstar(i) - pmns(i)) / pmns(i) wpls(i) = cpls(i) * sqrt (s1) wmns(i) = cmns(i) * sqrt (s2) enddo do i = 5, np5 zp = 4.d0 * taupls(i) * wpls(i) * wpls(i) zp = -zp * wpls(i) / (zp - gamp1 * (pstar(i) - ppls(i))) zm = 4.d0 * taumns(i) * wmns(i) * wmns(i) zm = zm * wmns(i) / (zm - gamp1 * (pstar(i) - pmns(i))) ustarp = upls(i) - (pstar(i) - ppls(i)) / wpls(i) ustarm = umns(i) + (pstar(i) - pmns(i)) / wmns(i) pstr1(i) = pstar(i) pstar(i) = pstar(i) 1 + (ustarm - ustarp) * (zm * zp) / (zm - zp) pstar(i) = max (smallp, pstar(i)) enddo c..test for convergence emax = 0.d0 do i = 5, np5 error = dabs (pstar(i) - pstr1(i)) / pstar(i) if (error .gt. emax) then emax = error imax = i endif enddo if (emax .le. tol) go to 60 enddo c..end of iteration loop 60 continue if (emax. gt. tol) then write (6,*) '***Warning*** Convergence failure in rieman.' write (6,*) 'Maximum error = ', emax, ' in zone ', imax endif c..calculate velocity do i = 5, np5 ustarp = upls(i) - (pstar(i) - ppls(i)) / wpls(i) ustarm = umns(i) + (pstar(i) - pmns(i)) / wmns(i) ustar(i) = 0.5d0 * (ustarp + ustarm) enddo return end subroutine tstep implicit none save c..compute new time step value using cfl condition c..bring in the common blocks include 'ppm.dek' c..local variables integer i double precision cflmax,s1,olddt,ceul cflmax = 0.d00 do i = 1, nr ceul = sqrt (gamma * press(i) / densty(i)) s1 = ceul / (rr(i) - rl(i)) if (s1 .lt. cflmax) go to 10 cflmax = s1 10 continue enddo olddt = dt dt = cfl / cflmax dt = min (dt, 1.2d0 * olddt) return end subroutine output implicit none save c..writes a result to a file c..bring in the common blocks include 'ppm.dek' c..local variables integer i double precision rreal(nd1) c..format statements 1000 format (1x,i6, 1p3e14.6, i6) 1005 format (i4, 5(2x, 1pe11.4)) open (3, file = 'results', status = 'unknown') write(3,1000) nstep, t, dt, gamma, nr do i = 1, nr write(3,1005) i, r(i), densty(i), vr(i), press(i), energy(i) enddo close(3) return end