C C define(zero,ifelse(ifdouble,1,0.d+00,0.0e+00)) define(one,ifelse(ifdouble,1,1.d+00,1.0e+00)) define(two,ifelse(ifdouble,1,2.d+00,2.0e+00)) define(three,ifelse(ifdouble,1,3.d+00,3.0e+00)) define(four,ifelse(ifdouble,1,4.d+00,4.0e+00)) define(five,ifelse(ifdouble,1,5.d+00,5.0e+00)) define(six,ifelse(ifdouble,1,6.d+00,6.0e+00)) define(seven,ifelse(ifdouble,1,7.d+00,7.0e+00)) define(nine,ifelse(ifdouble,1,9.d+00,9.0e+00)) define(ten,ifelse(ifdouble,1,10.d+00,1.0e+01)) define(twelve,ifelse(ifdouble,1,12.d+00,1.2e+01)) define(ffteen,ifelse(ifdouble,1,15.d+00,1.5e+01)) define(onept1,ifelse(ifdouble,1,1.1d+00,1.1e+00)) define(point9,ifelse(ifdouble,1,.9d+00,9.0e-01)) define(half,ifelse(ifdouble,1,.5d+00,5.0e-01)) define(qtr,ifelse(ifdouble,1,.25d+00,2.5e-01)) define(eighth,ifelse(ifdouble,1,.125d+00,1.25e-01)) define(fifth,ifelse(ifdouble,1,.2d+00,2.0e-01)) define(tenth,ifelse(ifdouble,1,.1d+00,1.0e-01)) define(sxy4th,ifelse(ifdouble,1,.015625d+00,1.5625e-02)) c define(klein,ifelse(ifdouble,1,1.0d-08,1.0e-08)) c subroutine ifelse(ifdouble,1,d_,) do_ppm98_deul_flux ^ ( xl, dx, ^ dalfac, darfac, facda, fcdazl, facdal, fa6dal, ^ fa6dar, flldal, flldar, frrdal, frrdar, ^ rho, p, u, ceul, p00, gamma1, unsmad, unsmup, ^ stpadb, rhoavl, ^ pavl, uavl, dvoll, dmassl, dmoml, detot1l, ^ courno, ssmlro, dt, N,nbdy_set) c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** ***** c ***** PPPPP ppppp m m ***** c ***** p p p p mm mm ***** c ***** ppppp ppppp m m m m ***** c ***** p p m m m ***** c ***** p p m m ***** c ***** p p m m ***** c ***** ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c c ************************************************************************** c * THIS CODE IS DO_PPM98_DEUL_FLUX. * c * THIS CODE WAS WRITTEN BY: WOODWARD RESEARCH GROUP. (c) 1998 * c * ALL RIGHTS RESERVED. VOID WERE PROHIBITED * c ************************************************************************** c c C WHAT This routine calculates the zone interface fluxes c for a direct eulerian hydrodynamic calculation in 1d. c c ASSUMPTIONS A Cartesian grid is assumed, which need not be uniform c 1D is assumed. (multi dimensions by directional splitting) c 1D is assumed. (multi dimensions by directional splitting) c A general metal) equation of state is used: C p(i) = p00(i) + gamma1(i) * rho(i) * ei(i) C C BOUNDARIES The necessary number of boundary zones to set (NBDY_SET) c is 4, three for interpolation, one for integration, c It will be further assumed that at the outset the c 2*NBDY_SET boundary zones [often referred to as fake c zones] contain valid values appropriate to the user's c chosen boundary conditions. This is the only means of c signalling boundary conditions to this routine The user c signals his understanding and acknowledges this by passing c to PPM98_DEUL0_FLUX the number of boundary zones set c on each side. This routine verifies that sufficient fake c zones have been provided before continuing. If not, an c advisory is displayed and the program terminated. C C DIMENSIONS c input It will be assumed that the dimensions of all zone-averaged c type quantities { e.g. RHO, P, U, RHONU, PNU, etc, ] are C [1-NBDY_SET:N_REAL_ZONES_IN+NBDY_SET]. C The dimensions of all [and I mean ALL] interface type C quantities [e.g. xl, dvoll, dmassl, pavl, uavl, anything C referenced to and interface, usually the lefthand zone c interface] are C [1-NBDY_SET:N_REAL_ZONES_IN+NBDY_SET+1]. C Note the "+1" please. There may be a quiz later. C C output At the completion of this entire routine, new zone-averaged C quantities [e.g. rhonu, pnu, unu] will have updated values C on the range C [1:N_REAL_ZONES_IN], C and any returned interface type thingies will have valid C values on the range [There are none at this time]: C [1:N_REAL_ZONES_IN+1]. c C c C DATA C input Input arrays are: C C dalfac: PPMLIB interpolation constant found from a previous c call to PPM98_INTERPOLATION_FACTORS C darfac: another PPMLIB interpolation constant C dacda: another PPMLIB interpolation constant C fcdazl: another PPMLIB interpolation constant C facdal: another PPMLIB interpolation constant C fa6dal: another PPMLIB interpolation constant C fa6dar: another PPMLIB interpolation constant C flldal: another PPMLIB interpolation constant C flldar: another PPMLIB interpolation constant C frrdal: another PPMLIB interpolation constant C frrdar: another PPMLIB interpolation constant C xl: the locations of the lefthand zone interfaces, c defining the computational grid. Please realize c that a numerical grid of n zones would require n+1 xl c values (interface locations.) C dx: The zone widths dx(i) = xl(i+1)-xl(i) C rho: The zone-averaged densities C p: The zone-averaged pressures C u: The zone-averaged (longitudinal) velocities C ceul: The zone-averaged eulerian sound speeds C p00: A zone-averaged equation of state coefficient C gamma1: A zone-averaged equation of state coefficient c A general metal) equation of state is used: C p(i) = p00(i) + gamma1(i) * rho(i) * ei(i) c c C Input scalers are: c C ssmlro: a trivial value for density C the above value is generally about C 1.0e-06 for single precision and C 1.0d-08 for dourble precision. C dt: the timestep. C n: The number zones (NOT counting fake zones) to be c processed. Upon completion of the timestep, c this many zones will have valid entries in the c output arrays. C nbdy_set: Number of boundary zones on either end set. This c is a check to ensure sufficient values have been c set by the calling routine. The number of c boundary zones required on either side must be c at least equat to the required values: 4 C C output Output arrays are: C unsmup: The "unsmoothness" indicator for the velocity c and pressure distributions. C unsmad: The "unsmoothness" indicator for the density c (adiab) distributions. C stpadb: The "steepening" or discontinuity indicator c for the density (adiab) distributions. c dt: timestep. C rhoavl: The time-averaged density at the lefthand c zone interface found from Riemann problems C pavl: The time-averaged pressure at the lefthand c zone interface found from Riemann problems C uavl: The time-averaged (longitudinal) velocity at c the lefthand zone interface found from Riemann c problems C dvoll: The volume advected across the lefthand interface C dmassl: The mass advected across the lefthand interface C dmoml: The (longitudinal) momentum advect across the c lefthand interface zone interface found from c Riemann problems C detot1l: The total energy (internal + longitudinal KE) c advected across the lefthand interface C courno: The Courant number (stability). To be stable c the domains of dependence for each zone interface c must lay within one zone. (That is, PPM is c explicit. Since a zone interface has domains of c of dependence on both the left and right, the c appropriate range of COURNO to examine is c 0:N+1, if the "real" data run from 1:N. c C C Output scaler is: C C There are none at this time. C C C c**************************************************************************** c**************************************************************************** c**************************************************************************** c**************************************************************************** c**************************************************************************** c**************************************************************************** c c ifelse(ifdouble,1,implicit real*8 (a-h,o-z), ) dimension xl(1-nbdy_set:N+nbdy_set+1), ^ dalfac(1-nbdy_set:N+nbdy_set), ^ darfac(1-nbdy_set:N+nbdy_set), ^ facda(1-nbdy_set:N+nbdy_set), ^ fcdazl(1-nbdy_set:N+nbdy_set), ^ facdal(1-nbdy_set:N+nbdy_set), ^ fa6dal(1-nbdy_set:N+nbdy_set), ^ fa6dar(1-nbdy_set:N+nbdy_set), ^ frrdal(1-nbdy_set:N+nbdy_set), ^ frrdar(1-nbdy_set:N+nbdy_set), ^ flldal(1-nbdy_set:N+nbdy_set), ^ flldar(1-nbdy_set:N+nbdy_set), ^ p00(1-nbdy_set:N+nbdy_set), ^ gamma1(1-nbdy_set:N+nbdy_set), ^ rho(1-nbdy_set:N+nbdy_set), ^ p(1-nbdy_set:N+nbdy_set), ^ u(1-nbdy_set:N+nbdy_set), ^ ceul(1-nbdy_set:N+nbdy_set), ^ unsmad(1-nbdy_set:N+nbdy_set), ^ unsmup(1-nbdy_set:N+nbdy_set), ^ stpadb(1-nbdy_set:N+nbdy_set), ^ courno(1-nbdy_set:N+nbdy_set), ^ uavl(1-nbdy_set:N+nbdy_set+1), ^ pavl(1-nbdy_set:N+nbdy_set+1), ^ rhoavl(1-nbdy_set:N+nbdy_set+1), ^ detot1l(1-nbdy_set:N+nbdy_set+1) dimension dx(1-nbdy_set:N+nbdy_set), ^ dvol(1-nbdy_set:N+nbdy_set), ^ dvoll(1-nbdy_set:N+nbdy_set+1), ^ dmassl(1-nbdy_set:N+nbdy_set+1), ^ dmoml(1-nbdy_set:N+nbdy_set+1) c c c c c c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** THIS CODE WAS WRITTEN BY ***** c ***** ***** c ***** PAUL R. WOODWARD ***** c ***** ***** c ***** May, 1992 ***** c ***** REVISED August, 1992 ***** c ***** REVISED September, 1992 ***** c ***** REVISED July, 1993 ***** c ***** REVISED June, 1997 ***** c ***** REVISED August, 1997 ***** c ***** BKEISED September, 1997 ***** c ***** ***** c ***** ***** c ***** ALL RIGHTS RESERVED ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c c c c c This routine computes the time-averaged density, pressure, and c velocity at each zone interface appropriate fora PPM Eulerian c hydrodynamic step. The time-averages are found using a Riemann c solver which treats all waves as discontinuities. c c c It is assumed that all boundary conditions have been implemented in c the calling routine by the proper addition of fake zones to both c ends of the input arrays. This code is 1d. c c Monotonicity constraints are applied to the hydrodynamical quantities only c where the functions are not somooth. c c c c c In June, 1993, this code was revised to improve its performance c on scalar microprocessors. This revision followed the example c of the code SUPERSCALARSPPM (6/14/93), in phys5062-93/week6. c For each interpolated quantity, measures of the lack of smooth- c ness of the function are constructed and used to govern the c application of monotonicity constraints, discontinuity steepen- c ing, and the computation of additional numerical dissipation. c c In August, 1997, this code was revised and extensively tested c so that it is ready for distribution to collaborators. c c In Februar 1998, many of these revisions were removed to that c normal people could read it. c allocatable :: rhol(:), pl(:), ul(:), 1 rhor(:), pr(:), ur(:), 2 drho(:), dp(:), du(:), 3 rho6(:), p6(:), u6(:) allocatable :: ullrp(:), pllrp(:), urlrm(:), prlrm(:), clagl(:) allocatable :: niterl(:) c n_required = 4 if ( nbdy_set .lt. n_required ) then call ppm98_message( ' ppm98_deul_flux_gamma Error. ***',0) call ppm98_message( ^ ' Note. the number of boundary zones set at either end of the ',0) call ppm98_message( ^ ' computational grid do not match the required number.',0) write(msg,6667) nbdy_set,n_required call ppm98_message( msg,0) 6667 format(' The calling program passed',i6,' boundary zones but ',i3, ' zones') write(msg,6668) N call ppm98_message( msg,-1) 6668 format(' are required. There were ',i9,' non-boundary zones passed.') endif c third = one / three forthd = four * third twelth = qtr * third c small = ifelse(ifdouble,1,1.0e-08,1.0d-08) smlfrc = ifelse(ifdouble,1,5.0d-03,5.0e-03) smlrho = ssmlro c allocate ( niterl(1-nbdy_set:N+nbdy_set+1), ^ rhol(1-nbdy_set:N+nbdy_set+1), ^ pl(1-nbdy_set:N+nbdy_set+1), ^ ul(1-nbdy_set:N+nbdy_set+1), 1 rhor(1-nbdy_set:N+nbdy_set+1), ^ pr(1-nbdy_set:N+nbdy_set+1), ^ ur(1-nbdy_set:N+nbdy_set+1), 2 drho(1-nbdy_set:N+nbdy_set), ^ dp(1-nbdy_set:N+nbdy_set), ^ du(1-nbdy_set:N+nbdy_set), 3 rho6(1-nbdy_set:N+nbdy_set), ^ p6(1-nbdy_set:N+nbdy_set), ^ u6(1-nbdy_set:N+nbdy_set)) allocate ( ullrp(1-nbdy_set:N+nbdy_set+1), ^ pllrp(1-nbdy_set:N+nbdy_set+1), ^ urlrm(1-nbdy_set:N+nbdy_set+1), 1 prlrm(1-nbdy_set:N+nbdy_set+1), 1 clagl(1-nbdy_set:N+nbdy_set+1)) KSTART = 1-nbdy_set KEND = N+nbdy_set KDO = N + 2*nbdy_set c c c Set up essential quantities. c CCC c do 1000 i = kstart,kend cesq = ceul(i)*ceul(I) c spdpls(i) = u(i) + ceul(i) c spdmns(i) = u(i) - ceul(i) c uabs = abs (u(i)) 1000 continue c c c Now find the spatially-averaged values in the left and right c domain of dependence for the lefthand zone interface for each zone. c c PPM98_EULER_DOMAIN_AVG calculates the left and right states for c a Riemann problem to each lefthand zone interface in the c following steps: c 1) Interpolate, with monotonicity constraints the c pressure and longitudinal velocity by constructing c interpolations of the +- Riemann invarients. c (Returns UL, UR, DU, U6, PL, PR, DP, P6, UNSMUP) c 2) Interpolate, with monotonicity constraints and c discontinuity detection, the density by c constructing interpolations of the 0 Riemann c invarient (streamline, or adiabatic invarient.) c (Returns UNSMAD, STPADB, RHOL, RHOR, DRHO, RHO6) c 3) Using the piecewise parabolic interpolations for c velocity, and pressure, find the averages in the c appropriate domains of dependence c c c c Now find the averages in the left and right c domain of dependence for the lefthand zone c interface for each zone. c c This routine, interpolates rho,u,p, finds the c averages of the appropriate dod, and returns these c averages: ullrp, pllrp, urlrm, prlrm. c c The average of density over the r0 DOD cannot c be found until the Riemann problem has been c solved. Therefore, the density interplations c are returned instead. nnbdy_set = 4 K = KSTART N = KDO - 2*nnbdy_set call ifelse(ifdouble,1,d_,) ppm98_Euler_domain_avg ^ ( dx(K), ^ dalfac(K), darfac(K), facda(K), fcdazl(K), facdal(K), fa6dal(K), ^ fa6dar(K), flldal(K), flldar(K), frrdal(K), frrdar(K), ^ rho(K), p(K), u(K), ceul(K), ^ ullrp(K), pllrp(K), urlrm(K), prlrm(K), clagl(K), ^ rhol(K), rhor(K), drho(K), rho6(K), ^ ul(K), ur(K), du(K), u6(K), ^ pl(K), pr(K), dp(K), p6(K), ^ unsmad(K), stpadb(K), unsmup(K), ^ courno(K), dt, N,nnbdy_set) kstart2 = kstart+nnbdy_set kend2 = kend-nnbdy_set kdo2 = kdo - 2*nnbdy_set c If one wished to examine the domain averages one could : c c print *,' i, ullrp(i), pllrp(i), urlrm(i), prlrm(i)' c do i=kstart2,kend2+1 c write(6,66) i, ullrp(i), pllrp(i), urlrm(i), prlrm(i) c enddo c c print *,' i, rhol(i), rho(i), rhor(i), drho(i), rho6(i)' c do i=kstart2,kend2+1 c write(6,66) i, rhol(i), rho(i), rhor(i), drho(i), rho6(i) c enddo c c c c *********************************************************************** c * * c * at this point, modifications to pllrp, prlrm to include body * c * forces, or to ullrp, urlrm to include geometrical effects * c * may be done. (See Colella and Woodward 1984). By modifying * c * the domain averages in this fashon, we may use a single * c * Riemann solver. * c * * c *********************************************************************** c c c Solve the Riemann problems and get the interface fluxes. c c c c PPM98_E_RIEMANN uses the domain averages of pressure and c velocity to solve Riemann's shocktube problem at each zone c interface, to provide time-averaged values of pressure c and velocity (returns PAVL and UAVL) at the zone interfaces. c The density pressure interpolations are then used to find c the time-averaged density (returns RHOAVL). From these c values, advected fluxes of volume, mass, momentum, and total c energy are computed. (Returns DVOLL, DMASSL, DMOML, and c DETOT1L) c nnbdy = 1 c KSTART2, KEND2, KDO2 refer to the dod interface values c c we require nnbdy additional zone on the left and right c of the density and pressure interpolations. We will use c these as our points of reference. Hence we adjust these \c c K values KSTART_RIEMANN = KSTART2-nnbdy KEND_RIEMANN = KEND2+nnbdy KDO_RIEMANN = KDO2+2*nnbdy c c NITERL is an integer array specifying the number of non-linear iterations c to take at each lefthand zone interface location. In smooth flow, or if c one intends to rely upon the later addition of additional dissapation, no c iterations are necessary. (NITERL(I)=0 for all i). However, situations c may arise inwich the user may wish to take additional iterations. Generally c only in the most extreme situations would more than 1 additional iteration c be required. The results from shock detection may be used to identify c such zones. Here, we rely upon additional dissapation. c K=KSTART_RIEMANN KEND_RIEMANN = KEND2+nnbdy N = KDO_RIEMANN - 2*nnbdy do i=k,k+n+1 niterl(i) = 0 enddo call ifelse(ifdouble,1,d_,) PPM98_E_RIEMANN ^ ( dx(K), ullrp(K), pllrp(K), urlrm(K), prlrm(K), ^ clagl(K), rhol(K), rho(K), rhor(K), drho(K), ^ rho6(K), pl(K), p(K), pr(K), dp(K), ^ p6(K), p00(K), gamma1(K), rhoavl(K), pavl(K), ^ uavl(K), dvoll(K), dmassl(K), dmoml(K), detot1l(K), ^ smlrho, dt, niterl(k), n, nnbdy) c c Update the zone counts c kstart3 = K + nnbdy kend3 = kend2 - nnbdy c c If one wished to examine the fluxes one could: c (Note, these are interface values, hence the +1) c print *,' i, rhoavl(i), pavl(i), uavl(i)' c do i=kstart3,kend3+1 c write(6,66) i, rhoavl(i), pavl(i), uavl(i) c enddo c c print *,' i, dvoll(i), dmassl(i), dmoml(i), detot1l(i)' c do i=kstart3,kend3+1 c write(6,66) i, dvoll(i), dmassl(i), dmoml(i), detot1l(i) c enddo 66 format(1h ,i3,1p,9(1x,e12.5)) c *************************************************************************** * * * c * The conservation laws may now be applied to obtain new zone-averaged * c * values of density, energy, and velocity. At this point, new pressures * c * are not possible because gases from adjacent zones have mixed. The * c * new equation of stae is not known. * c * * c *************************************************************************** c c c c c PPM call DONE c deallocate ( rhol, pl, ul, 1 rhor, pr, ur, 2 drho, dp, du, 3 rho6, p6, u6) deallocate ( niterl, ullrp, pllrp, urlrm, prlrm, clagl) return end