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_ppmde_1dc_paq_gamma ^ ( xl, rho, p, u, paq, rhonu, pnu, unu, ^ paqnu, eosgam, dt, smlrho, smallp, smallu, smalle, courmx, ^ nbyvol, nbymass, n_real_zones_in, nbdy_set) c c ************************************************************************** c * THIS CODE IS DO_PPMDE_PAQ_1DGAMMA. * 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 performs a 1D direct Eulerian hydrodynamic c update, with artificial diffusion. Provision is made c for passively advected quantities (PAQs). Both volume c weighted and mass weighted PAQs are supported. c c ASSUMPTIONS A Cartesian grid is assumed, which need not be uniform. c A gamma-law equation of state is used: C p(i) = ( eosgam - 1.0) * rho(I) * e_internal(i) C where eosgam [gamma] is the ratio of specific C heats. c Whatever effects the PAQs may have upon the physics c of the problem, it is assumed that the coupling c with the hydrodynamics is weak, so that the PAQs c may be treated as purely passive. C C C BOUNDARIES The necessary number of boundary zones to set (NBDY_SET) c is 5, three for interpolation, one for integration, and c Our Fiend Diffuse takes one. It will be further assumed c that at the outset the 2*NBDY_SET boundary zones [often c referred to as fake zones] contain valid values c appropriate to the user's chosen boundary conditions. C This is the only means of signalling boundary conditions C to this routine The user signals his understanding and c acknowledges this by passing to Do_PPM_DE_1DGAMMA c the number of boundary zones set on each side. This c routine verifies that sufficient fake zones have been c provided before continuing. If not, an advisory is c 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 c The PAQ arrays are two dimensional arrays which may be c considered arrays of arrays of zone-averaged quantities. c Their dimensions are [NXL:NXR,1:NA] where c NXL = 1-NBDY_SET c NXR = N_REAL_ZONES_IN+NBDY_SET, and c NA = NBYVOL+NBYMAS. c See Special note blow. c 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 For each PAQ ipaq, the valid range is C [1:N_REAL_ZONES_IN,ipaq]. 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 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 rho: The zone-averaged densities C p: The zone-averaged pressures C u: The zone-averaged velocities C paq: The zone-averaged Passively Advected Quantities. This c is assumed to be a TWO-DIMENSIONAL array. See special c note below. c C Input scalers are: C C eosgam: the ratio of specific heats (gamma) for a c gamma-law gas. C dt: the timestep. C smlrho: a trivial value for density C smallp: a trivial value for pressure C smallu: a trivial value for velocity C smalle: a trivial value for nergy C the above 4 values are generally about C 1.0e-06 for single precision and C 1.0d-08 for dourble precision. C nbyvol: The number of volume weighted PAQs. C nbymass: The number of mass weighted PAQs. C n_real_zones_in: 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: 5 C C output Output arrays are: C rhonu: The new zone-averaged densities at the new c timestep. C pnu: The new zone-averaged pressures at the new c timestep. C unu: The new zone-averaged velocities at the new c timestep. C paqnu: The new zone-averaged Passively Advected c Quantities. This is assumed to be a c TWO-DIMENSIONAL array. See special note below. C C C Output scaler is: C C courmx: The maximum Courant number found by this call C to this subroutine. C C C C SPECIAL NOTES c c paq The PAQ arrays paq and paqnu are two-dimensional arrays c of zone averaged quantities. The fastest running dimension c is the grid dimension, the second dimension is the paq c dimension. Thus one may think of them as arrays of arrays c of zone-averaged data. Specifically, the assumed dimensions c are c [1-NBDY_SET:N_REAL_ZONES_IN+NBDY_SET, 1:(NBYVOL+NBYMASS)]. c This assumption must be aheared to, or undesirable behaviour c may result. c c C VORSICHT HABEN! C C IT IS ASSUMED THAT APPROPRIATE VALUES OF THE DEPENDENT QUANTITIES C (xl, rho, p ,u) HAVE BEEN PLACED INTO THESE BOUNDARY ZONES C BY THE CALLING ROUTINE TO IMPLIMENT WHATEVER BOUNDARY CONDITIONS C THE USER DESIRED. C C VORSICHT HABEN! C C ALL, REPEAT _*ALL*_ (n_real_zones_in + 2*nbdy_set) C ZONES OF RHO, P, U, and (n_real_zones_in + 2**nbdy_set+1) C ZONES OF XL (note the +1) WILL BE USED AND MUST CONTAIN C VALID VALUES. C C VORSICHT HABEN! C C The user signals his understanding and acknowledges this by C passing to Do_PPM_DE_1DGAMMA the number of boundary zones C set on each side. Ve haf vays of checking dis. 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_real_zones_in + nbdy_set+1), ^ rho(1-nbdy_set:n_real_zones_in + nbdy_set), ^ p(1-nbdy_set:n_real_zones_in + nbdy_set), ^ u(1-nbdy_set:n_real_zones_in + nbdy_set), ^ rhonu(1-nbdy_set:n_real_zones_in + nbdy_set), ^ pnu(1-nbdy_set:n_real_zones_in + nbdy_set), ^ unu(1-nbdy_set:n_real_zones_in + nbdy_set) dimension paq(1-nbdy_set:n_real_zones_in + nbdy_set,nbyvol+nbymass), ^ paqnu(1-nbdy_set:n_real_zones_in + nbdy_set,nbyvol+nbymass) dimension ^ niterl(1-nbdy_set:n_real_zones_in + nbdy_set), ^ clagl(1-nbdy_set:n_real_zones_in + nbdy_set), ^ rhol(1-nbdy_set:n_real_zones_in + nbdy_set), ^ rhor(1-nbdy_set:n_real_zones_in + nbdy_set), ^ drho(1-nbdy_set:n_real_zones_in + nbdy_set), ^ rho6(1-nbdy_set:n_real_zones_in + nbdy_set), ^ pl(1-nbdy_set:n_real_zones_in + nbdy_set), ^ pr(1-nbdy_set:n_real_zones_in + nbdy_set), ^ dp(1-nbdy_set:n_real_zones_in + nbdy_set), ^ p6(1-nbdy_set:n_real_zones_in + nbdy_set), ^ ul(1-nbdy_set:n_real_zones_in + nbdy_set), ^ ur(1-nbdy_set:n_real_zones_in + nbdy_set), ^ du(1-nbdy_set:n_real_zones_in + nbdy_set), ^ u6(1-nbdy_set:n_real_zones_in + nbdy_set), ^ ullrp(1-nbdy_set:n_real_zones_in + nbdy_set), ^ pllrp(1-nbdy_set:n_real_zones_in + nbdy_set), ^ urlrm(1-nbdy_set:n_real_zones_in + nbdy_set), ^ prlrm(1-nbdy_set:n_real_zones_in + nbdy_set), ^ ceul(1-nbdy_set:n_real_zones_in + nbdy_set), ^ courno(1-nbdy_set:n_real_zones_in + nbdy_set), ^ detotl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ difuse_vel(1-nbdy_set:n_real_zones_in + nbdy_set), ^ difusl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dm(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dmasll(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dmasrl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dmassl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dminv(1-nbdy_set:n_real_zones_in + nbdy_set), ^ dmnew(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dmnu(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dmoml(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dvol(1-nbdy_set:n_real_zones_in + nbdy_set), ^ dvolinv(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dvoll(1-nbdy_set:n_real_zones_in + nbdy_set), ^ dx(1-nbdy_set:n_real_zones_in + nbdy_set), ^ smallda(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ unsmuta(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ steepa(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dal(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dvolfl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dmasfl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dmnuinv(1-nbdy_set:n_real_zones_in + nbdy_set), ^ dalfac(1-nbdy_set:n_real_zones_in + nbdy_set), ^ darfac(1-nbdy_set:n_real_zones_in + nbdy_set), ^ facda(1-nbdy_set:n_real_zones_in + nbdy_set), ^ fcdazl(1-nbdy_set:n_real_zones_in + nbdy_set), ^ facdal(1-nbdy_set:n_real_zones_in + nbdy_set), ^ fa6dal(1-nbdy_set:n_real_zones_in + nbdy_set), ^ fa6dar(1-nbdy_set:n_real_zones_in + nbdy_set), ^ frrdal(1-nbdy_set:n_real_zones_in + nbdy_set), ^ frrdar(1-nbdy_set:n_real_zones_in + nbdy_set), ^ flldal(1-nbdy_set:n_real_zones_in + nbdy_set), ^ flldar(1-nbdy_set:n_real_zones_in + nbdy_set), ^ dxinv(1-nbdy_set:n_real_zones_in + nbdy_set), ^ e_total(1-nbdy_set:n_real_zones_in + nbdy_set), ^ e_totnu(1-nbdy_set:n_real_zones_in + nbdy_set), ^ rhoavl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ stpadb(1-nbdy_set:n_real_zones_in + nbdy_set), ^ pavl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ uavl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ unsmad(1-nbdy_set:n_real_zones_in + nbdy_set), ^ unsmup(1-nbdy_set:n_real_zones_in + nbdy_set) character*100 msg CF90 allocatable :: courno(:), difuse_vel(:), CF90 ^ unsmad(:), unsmup(:), CF90 ^ stpadb(:), ceul(:), CF90 ^ rhoavl(:), pavl(:), CF90 ^ uavl(:), e_total(:), CF90 ^ e_totnu(:), dx(:), CF90 ^ dxinv(:), difusl(:), CF90 ^ dmassl(:), dmoml(:) CF90 allocatable :: dmasll(:), dmasrl(:), CF90 ^ detotl(:), dm(:), CF90 ^ dmnu(:), dmnew(:) CF90 allocatable :: dvol(:), dvoll(:) CF90 allocatable :: dvolinv(:), dminv(:) CF90 allocatable :: smallda(:), unsmuta(:), steepa(:), dal(:) CF90 allocatable :: dvolfl(:), dmasfl(:) CF90 allocatable :: dvolinv(:), dminv(:), dmnuinv(:) CF90 allocatable :: niterl(:), clagl(:), CF90 ^ rhol(:), pl(:), CF90 ^ rhor(:), pr(:), CF90 ^ rho6(:), p6(:), CF90 ^ drho(:), dp(:), CF90 ^ ul(:), ur(:), CF90 ^ du(:), u6(:), CF90 ^ ullrp(:), pllrp(:), CF90 ^ urlrm(:), prlrm(:) CF90 CF90 CF90 CF90 CF90 allocatable :: CF90 ^ dalfac(:), darfac(:), facda(:), fcdazl(:), CF90 ^ facdal(:), fa6dal(:), fa6dar(:), frrdal(:), CF90 ^ frrdar(:), flldal(:), flldar(:) CF90 small = klein gamma = eosgam c c We will create and use the following "K" indices as pointers to c the beginning and end of the working array. A count KDO will c also be maintained of the number of zones in the array to c process. Note that all of these shall refer to ZONE counts. c Quantities which depend upon lefthand zone interfaces, will c of course require one additional reference on the right. c the K*1's are appropriate whenever starting from scratch with c interpolation. KSTART1 = 1-NBDY_SET KEND1 = N_REAL_ZONES_IN + NBDY_SET KDO1 = N_REAL_ZONES_IN + 2 * NBDY_SET courmx = 0.0 n_required = 5 if ( nbdy_set .lt. n_required ) then call ppm98_message( ' Do_PPM_DE_PAQ_1DGAMMA 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) call ppm98_message( ^ ' This routine which performes a single step of a Direct Eulerian',0) call ppm98_message( ^ ' hydrodynamic update (including artificial diffusion) requires',0) write(msg,6666) n_required call ppm98_message( msg,0) 6666 format(1x, i3,' fake zones at each end of the array') call ppm98_message( ^ ' IT IS ASSUMED THAT APPROPRIATE VALUES OF THE DEPENDENT',0) call ppm98_message( ^ ' QUANTITIS (xl, rho, p ,u) HAVE BEEN PLACED INTO THESE',0) call ppm98_message( ^ ' BOUNDARY ZONES BY THE CALLING ROUTINE TO IMPLIMENT ',0) call ppm98_message( ^ ' WHATEVER BOUNDARY CONDITIONS THE USER DESIRED. ALL',0) call ppm98_message( ^ ' (n_real_zones_in + 2*nbdy_set)',0) call ppm98_message( ^ ' ZONES OF RHO, P, U and ',0) call ppm98_message( ^ ' (n_real_zones_in + 2*nbdy_set+1)',0) call ppm98_message( ^ ' ZONES OF XL (note the +1) WILL BE USED AND MUST ',0) call ppm98_message( ^ ' CONTAIN VALID VALUES.',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_real_zones_in call ppm98_message( msg,-1) 6668 format(' are required. There were ',i9,' non-boundary zones passed.') endif CF90 allocate ( CF90 ^ niterl(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ clagl(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ rhol(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ rhor(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ drho(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ rho6(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ pl(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ pr(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ dp(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ p6(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ ul(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ ur(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ du(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ u6(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ ullrp(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ pllrp(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ urlrm(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ prlrm(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ ceul(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ courno(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ detotl(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ difuse_vel(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ difusl(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dm(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dmasll(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dmasrl(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dmassl(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dminv(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ dmnew(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dmnu(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dmoml(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dvol(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ dvolinv(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dvoll(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ dx(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ smallda(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ unsmuta(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ steepa(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dal(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dvolfl(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dmasfl(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ dmnuinv(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 CF90 CF90 ^ dalfac(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ darfac(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ facda(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ fcdazl(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ facdal(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ fa6dal(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ fa6dar(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ frrdal(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ frrdar(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ flldal(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ flldar(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 CF90 ^ dxinv(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ e_total(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ e_totnu(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ rhoavl(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ stpadb(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ pavl(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ uavl(1-nbdy_set:n_real_zones_in + nbdy_set+1), CF90 ^ unsmad(1-nbdy_set:n_real_zones_in + nbdy_set), CF90 ^ unsmup(1-nbdy_set:n_real_zones_in + nbdy_set) CF90 ^ ) c c note that ALL interface thingies are dimensioned one more. c c c compute some essential quantities: c do i=kstart1, kend1 dx(i) = xl(I+1) - xl(i) dxinv(i) = 1.0 / dx(i) dvol(i) = dx(i) dvolinv(i) = dxinv(i) c NOTE these last two are true only in cartesian coords. dm(i) = rho(i) * dvol(i) dminv(i) = 1.0 / dm(i) ceul(i) = sqrt (eosgam * p(I) / rho(I) ) enddo c c Precompute various geometrical factors for use when interpolating c on a non-uniform grid. This permits greater efficiencies when c several quantities are to be interpolated. c c nbdy = 3 c nbdy_set should be <= NBDY_SET ! K = KSTART1 nzones = KDO1 - 2*nbdy c c precompute the interpolation factors call ifelse(ifdouble,1,d_,) ppm98_interpolation_factors ^ ( dx(K), dalfac(K), darfac(K), facda(K), fcdazl(K), ^ facdal(K), fa6dal(K), fa6dar(K), frrdal(K), frrdar(K), ^ flldal(K), flldar(K), nzones,nbdy) c Get the fluxes of the conserved hydrodynamical quantities: c These are mass, longitudinal momentum, and energy. c Get the fluxes of the conserved hydrodynamical quantities: c These are mass, longitudinal momentum, and energy. c This requires the following steps: c c 1) Interpolate, with monotonicity constraints the c pressure and longitudinal velocity by constructing c interpolations of the +- Riemann invarients. c (Returns 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) c 3) Averaging the piecewise parabolic interpolations for c density, velocity, and pressure, find the c averages in the appropriate domains of dependence c for a Riemann problem at each zone interface. c 4) Solve Riemann shock tube problems at each lefthand c zone interface to obtain appropriate time-averaged c values for density, pressure, and longitudinal c momentum at these interfaces. c (Returns RHOAVL, PAVL, UAVL.) c 5) From these time averages, construct interface fluxes c (in the Eulerian sense) of mass, momentum, and c total energy. c (Returns DMASSL, DMOML, and DETOTL.) c 6) Determine the Courant condition (stability) for c each _INTERFACE_. Simply put, in order for c the method to be stable, the various Domains c of dependence for the Riemann problems c cannont exceed the width of either the zone c to the left or the zone to the right. This c make the courant number array COURNO appear c to be an interface-like quantity. Its valid c range shall be [K:K+NZONES] c c c third = one / three forthd = four * third twelth = qtr * third smlfrc = ifelse(ifdouble,1,5.0d-03,5.0e-03) c N = n_real_zones_in KSTART = 1-nbdy_set KEND = N+nbdy_set KDO = N + 2*nbdy_set 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 Note: c The average of density over the r0 DOD cannot c be found until the Riemann problem has been c solved. Therefore, the density interplotions c are returned only. 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), 6 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) c c Update the zone counts c 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 print *,' i, rhol(i), rho(i), rhor(i), drho(i), rho6(i)' c do i=kstart2-1,kend2+1 c write(6,66) i, rhol(i), rho(i), rhor(i), drho(i), rho6(i) c enddo 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 PPM98_E_RIEMANN_GAMMA 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. c Generally only in the most extreme situations would more than 1 c additional iteration be required. The results from shock detection may c be used to identify such zones. Here, we rely upon additional c dissapation. c K=KSTART_RIEMANN N = KDO_RIEMANN - 2*nnbdy do i=k,k+n+1 niterl(i) = 0 enddo call ifelse(ifdouble,1,d_,) PPM98_E_RIEMANN_GAMMA ^ ( 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), rhoavl(K), pavl(K), uavl(K), dvoll(K), ^ dmassl(K), dmoml(K), detotl(K), gamma, smlrho, ^ smallp, dt, niterl(K), n, nnbdy) c c Update the zone counts c kstart3 = kstart2 kend3 = kend2 kdo3 = kdo2 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), detotl(i)' c do i=kstart3,kend3+1 c write(6,66) i, dvoll(i), dmassl(i), dmoml(i), detotl(i) c enddo c 66 format(1h ,i3,1p,9(1x,e14.7)) c c c *************************************************************************** * * * c * The conservation laws may now be applied to obtain new zone-averaged * c * values of density, pressure (or energy), and velocity. * c * * c *************************************************************************** c c c******************************************************************* c the K*3 will be appropriate for working with fluxes (before diffuse) c find the courant number from the above step. When a routine is c called with an array such as COURNO, the old values of COURNO are c overwritten. Thus the Courant number should be determined after each c such routine, or different arrays used. c c To be stable when a flux at xl(i), the left hand edge of zone I, is c computed, the courant numbers in zone I-1 and I need to be <= 1.0 c Thus we examine Courant numbers over the KSTART3-1,KEND3+1 c courmxh = courno(KSTART3 -1) do i=KSTART3 , KEND3+1 courmxh = max (courmxh, courno(i)) enddo c c in order to update the zones, which is done with *conserved quantities* c (mass, momenta, energy), we need to compute the old zone averaged total c energy, which was not provided as initial data. How rude. It is c computed only over the range required. c gamm1i = 1.0 / ( eosgam - 1.0) do i=KSTART2 ,KEND2 ek = 0.5 * u(I)*u(I) eint = gamm1i * p(i) / rho(i) e_total(i) = eint + ek enddo c c c Using the above fluxes, update the zone averages of conserved quantities c K = KSTART2 nzones = KDO2 call ifelse(ifdouble,1,d_,)ppm98_deul_update_cart ^ ( dx(K), dvol(K), dm(K), rho(K), u(K), ^ e_total(K), dmassl(K), dmoml(K), pavl(K), detotl(K), ^ rhonu(K), unu(K), e_totnu(K), dmnu(K), smlrho, ^ smallp, smallu, smalle, dt, nzones) c no fake zones necessary in the above routine, so the "K" count is c unchanged. The updated zone-averages are valid on the range c kstart2,kend2, and could be examined by: c print *,' Finished with the hydro step.' c print *,' i, rhonu(i), unu(i), e_totnu(i), pnu(i)' c do i=kstart2,kend2 c ekk = 0.5 * (unu(i)**2) c eii = e_totnu(i) - ekk c peenu = (eosgam-1.0) * rhonu(i)*eii c write(6,666) i, rhonu(i), unu(i), e_totnu(i), peenu c enddo c c Now we have the information necessary to start on the PAQs. c if (nbyvol .eq. 0) goto 5000 c c First construct volume weightings: c c NOTE, for this example I don't really know what A is, so I will use c SMALL for smallda. If A is a velocity (e.g. transverse c velocity) then small * (eulerian sound speed) would be appropriate. c A positive definite quantiity such as moisture content, then c smallda = small(I)*a(i) but if a can be zero (or negative, then c this latter choice would be inappropriate. Apparently, c smallda should be an input quantity, and would thus be another 2D c array. That is too bad. c c We will use smallda for each PAQ, thus only one 1d strip array is c required. c do i=KSTART1,KEND1 smallda(i) = small enddo c c Obtain advected volume fractions c do i=KSTART2,KEND2+1 if( dvoll(i) .le. 0.0 ) then dvolfl(i) = -dvoll(i) * dvolinv(i) else dvolfl(i) = dvoll(i) * dvolinv(i-1) endif dvolfl(i) = min (one, dvolfl(i)) enddo c Advecte the individual PAQS: do ipaq = 1,nbyvol c c K = KSTART1 c c nbdy_za is the number of fake zones of zone-averaged data nbdy_za = 4 c nbdy_zl is the number of fake zones of zone-interface data nbdy_zl = 1 NZONES = KDO1 -2*nbdy_za c PPM98_ADVECT_PASSIVE0 performs the entire advection step. Its c steps are: c 1) Interpolate, with discontinuity detection and c monotonicity constraints. (Returns UNSMUTA c and STEEPA) c 2) Integrate the piecewise parabolic interpolations c over the appropriate fraction of each zone. c 3) Construct weighted fluxes from step 2) using the c provided interface weightings, here dvoll, c the volume advected at the lefthand zone interface. c (Returns DAL) c 4) Update the zone averages by conservatively differencing c the interface fluxes from step 3) and using the c old zone volume (DVOL) and the inverse of the new c zone volume (DVOLINV) as weightings. (For this c Eulerian scheme the old and new zone volumes are c the same and DVOLINV(i) = 1.0/DVOL(i).) (Returns c PAQNU). c c These individual steps are also available as separate c modules. c call ifelse(ifdouble,1,d_,) ppm98_advect_passive ^ ( dalfac(K), darfac(K), facda(K), fcdazl(K), facdal(K), fa6dal(K), ^ fa6dar(K), flldal(K), flldar(K), frrdal(K), frrdar(K), ^ dvolfl(K), dvoll(K), dvol(K), ^ dvolinv(K), paq(K,ipaq), smallda(K), ^ unsmuta(K), steepa(K), dal(K), ^ paqnu(K,ipaq), nzones, nbdy_za, nbdy_zl ) c enddo 5000 continue c Continue with the mass weighted PAQs (if any) c First construct mass weightings: if (nbymass .eq. 0) goto 5100 c c c We will require inverses of both the old and new zone masses: c do i=KSTART2-1,KEND2+2 dminv(i) = 1.0/dm(i) enddo do i=KSTART2,KEND2 dmnuinv(i) = 1.0/dmnu(i) enddo c c see note above about smallda in the volume section c do i=KSTART2,KEND2+1 smallda(i) = small enddo c c Obtain advected mass fractions c do i=KSTART2,KEND2+1 if( dmassl(i) .le. 0.0 ) then dmasfl(i) = -dmassl(i) * dminv(i) else dmasfl(i) = dmassl(i) * dminv(i-1) endif dmasfl(i) = min (one, dmasfl(i)) enddo c Advecte the individual PAQS: do ipaq = nbyvol+1,nbyvol+nbymass c K = KSTART1 c c nbdy_za is the number of fake zones of zone-averaged data nbdy_za = 4 c nbdy_zl is the number of fake zones of zone-interface data nbdy_zl = 1 NZONES = KDO1 -2*nbdy_za c PPM98_ADVECT_PASSIVE performs the entire advection step. Its c steps are: c 1) Interpolate, with discontinuity detection and c monotonicity constraints. (Returns UNSMUTA c and STEEPA) c 2) Integrate the piecewise parabolic interpolations c over the appropriate fraction of each zone. c 3) Construct weighted fluxes from set 2) using the c provided interface weightings, here dmass, c the mass advected across the lefthand zone interface. c (Returns DAL) c 4) Update the zone averages by conservatively differencing c the interface fluxes from step 3) and using the c old zone mass (DM) and the inverse of the new c zone mass (DMNUINV) as weightings. (Recall that c the zone masses have changed during the c hydrodynamic step above.) c c These individual steps are also available as separate c modules. c c call ifelse(ifdouble,1,d_,) ppm98_advect_passive ^ ( dalfac(K), darfac(K), facda(K), fcdazl(K), facdal(K), fa6dal(K), ^ fa6dar(K), flldal(K), flldar(K), frrdal(K), frrdar(K), ^ dmasfl(K), dmassl(K), dm(K), ^ dmnuinv(K), paq(K,ipaq), smallda(K), ^ unsmuta(K), steepa(K), dal(K), ^ paqnu(K,ipaq), nzones, nbdy_za, nbdy_zl ) c enddo 5100 continue c C ********************************************************** C * * C * Done with the hydro step, begin diffusion step * C * * C ********************************************************** c c zone centerred difusion velocities are required on the range c zone zone just updated. Any bondary requirements (offsets) c are built in into the offsets in the calling statement. c c Therefore the correct relative beginning locations is KSTART2 c and the number of zone to process is KDO2 c c BUT REALIZE THAT PPM98_EUL_DIFFUSE_1D OPERATES ON THE INITIAL C DATA. c c KSTARTDF = KSTART2 KENDDF = KEND2 KDODF = KDO2 c c c Get difusion velocities K=KSTARTDF nzones = KDODF call ifelse(ifdouble,1,d_,) ppm98_eul_difuse_1d ^ ( dx(k-2), dx(k-1), dx(K), dx(k+1), dx(k+2), ^ u(k-2), u(k-1), u(K), u(k+1), u(k+2), p(k-2), p(K), p(k+2), ^ rho(k-2), rho(K), rho(k+2), ceul(k-2), ceul(K), ceul(k+2), ^ difuse_vel(K), smallu, nzones) c update the zone counts KSTARTDF1 = KSTARTDF KENDDF1 = KENDDF KDODF1 = KDODF c At this point, we should have one extra diffusion velocity c at each end of the working array. These are ZONE-CENTERED c velocities, from which INTERFACE centerred fluxes will c be constructed. This construction requires this extra fake zone. C Now apply diffusion c c STEP 1: Get the diffusive zone interface fluxes c from the zone centerred difusion velocities nbdy_difuse_fluxes = 1 K = KSTARTDF1 nzones = KDODF1 - 2*nbdy_difuse_fluxes call ifelse(ifdouble,1,d_,) ppm98_difuseflux_hydro_cart ^ ( dx(K), difuse_vel(K), ^ rhonu(K), unu(K), e_totnu(K), ^ dvoll(K), dmasll(K), dmasrl(K), ^ dmassl(K), dmoml(K), detotl(K), ^ courno(K), dt, nzones, ^ nbdy_difuse_fluxes ) KSTARTDF2 = KSTARTDF1 + nbdy_difuse_fluxes KENDDF2 = KENDDF1 - nbdy_difuse_fluxes courdf2 = courno(KSTARTDF2) do i=KSTARTDF2+1,KENDDF2 courdf2 = max (courdf2,courno(i)) enddo c STEP 2: apply the diffusive fluxes c c (Final update:, whew!) c c We occasionally get underflows by allowing absurdly small c velocities to be computed. Therefore, we will set such c velocities to zero below. c ssmalu = small * smallu ssmlu2 = ssmalu * ssmalu c gamma1 = eosgam - 1.0 do 9100 i = KSTARTDF2,KENDDF2 c c save old dmnu to cherish later: dmnew(i) = dmnu(i) c this saves a flop or two: if ((dvoll(i) + dvoll(i+1)) .ne. zero) then c dmnew(i) = rhonu(i) * dvol(I) dmnu(i) = dmnew(i) + dmassl(i) - dmassl(i+1) rhonu(i) = dmnu(i) * dvolinv(i) dmnui = one / dmnu(i) unu(i) = (unu(i) * dmnew(i) + dmoml(i) - dmoml(i+1)) * dmnui e_totnu(i) = (e_totnu(i) * dmnew(i) + ^ detotl(i) - detotl(i+1)) * dmnui 9099 continue endif usq = unu(i) * unu(i) c if ((unu(i) * unu(i)) .lt. ssmlu2) unu(i) = 0.0 if ( usq .lt. ssmlu2) then unu(i) = 0.0 usq = 0.0 endif c eknu = 0.5 * (unu(i) * unu(i)) eknu = 0.5 * usq e_intnu = max (smalle, (e_totnu(i) - eknu)) pnu(i) = gamma1 * rhonu(i) * e_intnu c write(6,666) i, rhonu(i), unu(i), e_totnu(i), pnu(i) 9100 continue c c Now difuse the passive quantities c c First do the volume weighted paqs if (nbyvol .eq. 0) goto 6000 K = KSTARTDF1 nbdy = 1 nzones = KDODF1 - 2*nbdy do ipaq = 1,nbyvol c c Use the diffusion volumes DVOLL at the lefthand zone interface to c find the volume weighted interface difusion fluxes DAL call ifelse(ifdouble,1,d_,) ppm98_difuseflux_passive ^ (paqnu(K,ipaq), dvoll(K), dvoll(K), dal(K), nzones, nbdy) KSTARTDF = KSTARTDF1 + nbdy KENDDF = KENDDF1 - nbdy c c Now update these quantities by conservative differencing the c volume weighted interfaces fluxes found above. c Could call ppm98_weighted_update here, but it is assumed that c most of the time dvoll(I) is zero, since difusion is applied c only in shocks. Thus this loop is a little more efficient. c when large number of PAQs are involved, it would be best to c construct a flag equalt to the sum of dvoll(i) + dvoll(i+1). c Clarity of proceedure dictates otherwise here. do 6800 i = KSTARTDF,KENDDF ghoti = paqnu(i,ipaq) if ((dvoll(i) + dvoll(i+1)) .ne. zero) then paqnu(i,ipaq) = paqnu(i,ipaq) + (dal(i) - dal(i+1)) * dvolinv(i) endif 6800 continue enddo 6000 continue c Now do the mass weighted paqs if (nbymass .eq. 0) goto 7000 K = KSTARTDF1 nbdy = 1 nzones = KDODF1 - 2*nbdy KSTARTDF = KSTARTDF1 + nbdy KENDDF = KENDDF1 - nbdy do ipaq = nbyvol+1,nbyvol+nbymass c c Use the leftward and rightward diffusion volumes (MASLL, MASRL) at the c lefthand zone interface to find the mass weighted interface difusion c fluxes DAL call ifelse(ifdouble,1,d_,) ppm98_difuseflux_passive ^ (paqnu(K,ipaq), dmasll(K), dmasrl(K), dal(K), nzones, nbdy) KSTARTDF = KSTARTDF1 + nbdy KENDDF = KENDDF1 - nbdy c c Now update these quantities by conservative differencing the c volume weighted interfaces fluxes found above. Again, c could call ppm98_weighted_update here, but it is assumed that c most of the time dvoll(I) is zero, since difusion is applied c only in shocks. Thus this loop is a little more efficient. c when large number of PAQs are involved, it would be best to c construct a flag equalt to the sum of dvoll(i) + dvoll(i+1). c Clarity of proceedure dictates otherwise here. do 7800 i = KSTARTDF,KENDDF ghoti = paqnu(i,ipaq) if ((dvoll(i) + dvoll(i+1)) .ne. zero) then dmnui = one / dmnu(i) c here we have used thecherished masses dmnew! paqnu(i,ipaq) = (paqnu(i,ipaq) * dmnew(i) + ^ dal(i) - dal(i+1)) * dmnui endif 7800 continue c enddo 7000 continue c c The final Courant number is the maximum of the c hydrodynamics Courant number and the diffusive c Courant number. c courmx = max (courmxh, courdf2 ) CF90 deallocate (courno, difuse_vel) CF90 deallocate ( unsmad, unsmup, stpadb, ceul) CF90 deallocate ( rhoavl, pavl, uavl) CF90 deallocate ( e_total, e_totnu, dx, dxinv) CF90 deallocate ( difusl, dmassl, dmoml, detotl) CF90 deallocate ( dmasll, dmasrl) CF90 CF90 deallocate (dvol, dvoll) CF90 deallocate (dm, dmnu, dmnew ) CF90c CF90 deallocate (smallda, unsmuta, steepa, dal , CF90 ^ dvolfl, dmasfl) CF90 deallocate (dvolinv, dminv, dmnuinv) CF90 CF90 deallocate (niterl, clagl) CF90c CF90 CF90 deallocate (dalfac, darfac, facda, fcdazl, CF90 ^ facdal, fa6dal, fa6dar, frrdal, CF90 ^ frrdar, flldal, flldar) CF90 CF90 deallocate ( niterl, CF90 ^ rhol, pl, CF90 ^ rhor, pr, CF90 ^ rho6, p6, CF90 ^ drho, dp, CF90 ^ ul, ur, CF90 ^ du, u6, CF90 ^ ullrp, pllrp, CF90 ^ urlrm, prlrm) return end