c subroutine ifelse(ifdouble,1,d_,)do_ppmlr0_2dc_gamma ^ ( xl, rhozbb, rho, rhoztt, pzbb, p, pztt, ^ uzbb, u, uztt, utzbb, utzb, ut, utzt, ^ utztt, rhonu, pnu, unu, utnu, eosgam, dt, ^ smlrho, smallp, smallu, smalle, courmx, n_real_zones_in, ^ nbdy_set) c c c ************************************************************************** c * THIS CODE IS DO_PPMLR0_2DC_GAMMA. * c * THIS CODE WAS WRITTEN BY: WOODWARD RESEARCH GROUP. (c) 1999 * c * ALL RIGHTS RESERVED. * c ************************************************************************** c C WHAT This routine performs a single pass of a 2D c hydrodynamic calculation in two parts: c a Lagrangian step followed by a remap back c to the original grid. c c ASSUMPTIONS A Cartesian grid is assumed, which is initially 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 C BOUNDARIES The necessary number of boundary zones to set (NBDY_SET) c is 9, three for interpolation, one for integration, in both c the Lagrangian step and remap step for a subtotal of 8, 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_LR0_2DC_GAMMA 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 abruptly 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 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. 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 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 rhozbb: The zone-averaged densities two rows below the pass c of interest. C rhoztt: The zone-averaged densities two rows above the pass c of interest. C p: The zone-averaged pressures in the pass of interest C pzbb: The zone-averaged pressures two rows below the pass c of interest. C pztt: The zone-averaged pressures two rows above the pass c of interest. C u: The zone-averaged longitudinal-velocities in the pass c of interest C uzbb: The zone-averaged longitudinal-velocities two rows c below the pass of interest. C uzb: The zone-averaged longitudinal-velocities one row c below the pass of interest. C uzt: The zone-averaged longitudinal-velocities one row c above the pass of interest. C uztt: The zone-averaged longitudinal-velocities two rows c above the pass of interest. C ut: The zone-averaged transverse-velocities in the pass of c interest C utzbb: The zone-averaged transverse-velocities two rows below c the pass of interest. C utzb: The zone-averaged transverse-velocities one row below c the pass of interest. C utzt: The zone-averaged transverse-velocities one row above c the pass of interest. C utztt: The zone-averaged transverse-velocities two rows above c the pass of interes 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 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 value: 9. 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 utnu: The new zone-averaged transverse- velocities c at the new timestep. 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 directional splitting c c PPM performs multi-dimensional hydrodynamics by directional c splitting, in which a one-dimensional operator is applied c in alteration, e. g. a sequence of passes x-y-y-x. During c the hydrodynamic step, the transverse velocities are c advected passively. Therefore the following proceedure c will be followed: c c 1) Find the the fluxes of the active hydrodynamical c quantities using ppm98_deul0_flux_gamma_1d. The c active hydro. quantities are: rho, pressure, c longitudinal velocity (u). c 2) using these fluxes to update the zone averages of c density (rhonu), longitudinal velocity (unu), and c a partial update of the total energy (e_totnu) c using ppm98_deul_update_1d. c c The portion of the total energy updated is the c internal energy + the longitudinal kinetic energy. c That is we will write the total energy as c c ET(I) = EI(I) + EKL(I) + EKT(I), c where c c ET(I) = total specific energy in zone I c c EI(I) = specific internal energy in zone I c c EKL(I) = total specific longitudinal kinetic c energy in zone I c ekl(i) = 0.5 * u(i)*u(i) c c EKT(I) = total specific transverse kinetic c energy in zone I c ekt(i) = 0.5 * ut(i)*ut(i). c c 3) From the fluxes at the zone interfaces combined c with the old and new zone masses, it is possible c to find the new zone averages of the transverse c velocity (UTNU). The fluxes of the transverse c velocity are then used to determine the transverse c kinetic energy portion of the total energy, which c is then added to the earlier partial ET update of c step2. c C C C C************************************************************************* C************************************************************************* C************************************************************************* C************************************************************************* C************************************************************************* C************************************************************************* C************************************************************************* C************************************************************************* C************************************************************************* C C these permit switching btw double and single precision 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************************************************************************* 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), ^ rhozbb(1-nbdy_set:n_real_zones_in + nbdy_set), ^ rho(1-nbdy_set:n_real_zones_in + nbdy_set), ^ rhoztt(1-nbdy_set:n_real_zones_in + nbdy_set), ^ pzbb(1-nbdy_set:n_real_zones_in + nbdy_set), ^ p(1-nbdy_set:n_real_zones_in + nbdy_set), ^ pztt(1-nbdy_set:n_real_zones_in + nbdy_set), ^ uzbb(1-nbdy_set:n_real_zones_in + nbdy_set), ^ u(1-nbdy_set:n_real_zones_in + nbdy_set), ^ uztt(1-nbdy_set:n_real_zones_in + nbdy_set), ^ utzbb(1-nbdy_set:n_real_zones_in + nbdy_set), ^ utzb(1-nbdy_set:n_real_zones_in + nbdy_set), ^ ut(1-nbdy_set:n_real_zones_in + nbdy_set), ^ utzt(1-nbdy_set:n_real_zones_in + nbdy_set), ^ utztt(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), ^ utnu(1-nbdy_set:n_real_zones_in + nbdy_set) character*100 msg allocatable :: denl(:) allocatable :: dvolfl(:) allocatable :: courno(:), difuse_vel(:), ^ unsmad(:), unsmup(:), ^ unsmutut(:), unsmooth(:), ^ stpadb(:), ceul(:), ^ rhoavl(:), pavl(:), ^ uavl(:), e_total(:), ^ e_totnu(:), dx(:), ^ dxinv(:), difusl(:), ^ dmassl(:), dmoml(:) allocatable :: xllag(:), rholag(:), ^ plag(:), ulag(:), ^ e_totlag(:), dxlag(:) allocatable :: dmasll(:), dmasrl(:), ^ detotl(:), dm(:), ^ dmnu(:), dmnew(:) allocatable :: dvol(:), dvoll(:) allocatable :: dvolinv(:), dminv(:) allocatable :: ^ smalldu(:), steeput(:), dmomtl(:), shockd(:), ^ dmassfl(:), dmnuinv(:), dektl(:), ^ dalfac(:), darfac(:), facda(:), fcdazl(:), ^ facdal(:), fa6dal(:), fa6dar(:), frrdal(:), ^ frrdar(:), flldal(:), flldar(:) small = klein 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 = 9 if ( nbdy_set .lt. n_required ) then call ppm98_message( ' Do_PPM_LR0_2DC_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) call ppm98_message( ^ ' This routine which performes a single step of a Lagrangian+Remap',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 call ppm98_message( msg,-1) 6668 format(' are required. There were ',i9,' non-boundary zones passed.') stop endif allocate ( ^ ceul(1-nbdy_set:n_real_zones_in + nbdy_set), ^ courno(1-nbdy_set:n_real_zones_in + nbdy_set), ^ denl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dektl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dmnuinv(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), ^ dmassfl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dvolfl(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), ^ smalldu(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ steeput(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dmomtl(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ shockd(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dx(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), ^ unsmutut(1-nbdy_set:n_real_zones_in + nbdy_set), ^ unsmooth(1-nbdy_set:n_real_zones_in + nbdy_set), ^ unsmad(1-nbdy_set:n_real_zones_in + nbdy_set), ^ unsmup(1-nbdy_set:n_real_zones_in + nbdy_set) ^ ) allocate ( ^ xllag(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ dxlag(1-nbdy_set:n_real_zones_in + nbdy_set+1), ^ rholag(1-nbdy_set:n_real_zones_in + nbdy_set), ^ plag(1-nbdy_set:n_real_zones_in + nbdy_set), ^ ulag(1-nbdy_set:n_real_zones_in + nbdy_set), ^ e_totlag(1-nbdy_set:n_real_zones_in + nbdy_set)) c c note that ALL interface thingies are dimensioned one more. 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 K = KSTART1 nzones = KDO1 - 2*nbdy c c Get the fluxes of the conserved hydrodynamical quantities: c These are mass, longitudinal momentum, and energy. Note, c the transverse velocity is merly passively advected. c do_ppm98_lag_flux_gamma requires 4 fake zones on each end of grid c c c c DO_PPM98_LR_FLUX_GAMMA produces the Lagrangian zone interface fluxes in the c longitudinal direction (the direction of this pass) for the c conserved hydrodynamic quantities mass, longitudinal momentum c and total energy in 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) Using these piecewise parabolic interpolations for c , velocity and pressure, find the c averages in the appropriate domains of dependence c for a Riemann problem at each zone interface. c 3) Solve Riemann shock tube problems at each lefthand c zone interface to obtain appropriate time-averaged c values for pressure, and longitudinal velocity at c these interfaces. Time averages of density are not c required for a Lagrangian calculation. c (Returns PAVL, UAVL.) c 4) 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 These individual steps are also available as separate c modules. c nbdy = 4 K = KSTART1 NZONES = KDO1 -2*nbdy c NZONES is the count of the number of zones which will have valid c entries in the output from DO_PPM98_LAG0_FLUX_GAMMA c call ifelse(ifdouble,1,d_,)do_ppm98_lag0_flux_gamma ^ ( xl(K), dx(K), ^ rho(K), p(K), u(K), ceul(K), unsmup(K), ^ pavl(K), uavl(K), courno(K), eosgam, smallp, ^ dt, nzones, nbdy) c update the zone counts. These are the valid counts for the interface c fluxes, dmassl, dvoll, etc. The number of these boundary fluxes which c are for zone THIS program considers fake are (beginning #) - (used #) c This program began with nbdy_set, and nbdy were just used. KSTART2 = KSTART1 + nbdy KEND2 = KEND1 - nbdy KDO2 = KDO1 - 2*nbdy c c Thus the fluxes are now valid on the range (kstart2,kend2+1. To examine c them one could do: c print *,' i, pavl(i), uavl(i)' c do i=kstart2,kend2+1 c print *,' i=',i c write(6,666) i, pavl(i), uavl(i) c enddo 666 format(1h ,i4,1p,9(1x,e10.3)) 66 format(1h ,i4,1p,9(1x,e10.3)) c the K*2 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 KSTART2-1,KEND2+1 c courmxh = courno(KSTART2 -1) do i=KSTART2 , KEND2+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. Because of the c directional splitting, we will use only the internal + longitudinal c portion of the kinetic energy. The transverse portion of the kinetic c energy will be updated when we have the transverse momentum flux. This c final term will then be added to the new zone averaged energy found here. c c gamm1i = 1.0 / ( eosgam - 1.0) do i=KSTART2 ,KEND2 ekl = 0.5 * u(I)*u(I) eint = gamm1i * p(i) / rho(i) e_total(i) = eint + ekl enddo c c c PPM98_LAG_UPDATE_CART upates the zone averages by conserviatvely c differencing the fluxes found above, assuming cartesian c coordinates and no body forces. Body forces would c entail and additional term added to the new zone-averages of c velocity and total energy AFTER this routine is called. c c The suffix "lag" will be appended to the updated quantities C (XLLAG, DXLAG, RHOLAG, ULAG, E_TITLAG) to indicate that these c are the quantities after the lagrangian step. (DMLAG == DM of c course) Traditional, quantities were renamed at this point, c over-writing the initial data. This is now avoided, as well as c the gratuitous copying. c K = KSTART2 nzones = KDO2 call ifelse(ifdouble,1,d_,)ppm98_lag_update_cart ^ ( xl(K), dx(K), dm(K), u(K), ^ e_total(K), uavl(K), pavl(K), xllag(k), dxlag(K), ^ rholag(K), ulag(K), 6 e_totlag(K), smallu, smalle, dt, nzones) c c 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 c c c print *,' Finished with the lag step.' c print *,' i, xllag(i), rholag(i), ulag(i), e_totlag(i), plag(i)' do i=kstart2,kend2 ekk = 0.5 * (ulag(i)**2) eii = e_totlag(i) - ekk plag(i) = (eosgam-1.0) * rholag(i)*eii c write(6,666) i,xllag(i), rholag(i), ulag(i), e_totlag(i), plag(i) enddo c C ********************************************************** C * * C * Done with the Lagrangian step, begin remap step * C * * C ********************************************************** C c The grid is now definitely non uniform. So, first pre-compute c the geometrical factors (the old factors, if any, are destroyed) c c the interpolation routine need 3 zones: nbdy = 3 c c in the interpolation routines, we require a value for the zone-averaged c quantity being interpolated which can be considered negligable. This c is used to redece the effects of round off and numerical noise. For c velocity, an appropriate value is small*(eulerian sound speed) do i=KSTART2 ,KEND2 smalldu(i) = small * ceul(i) enddo K = kstart2 nzones = KDO2 - 2*nbdy c c precompute the interpolation factors call ifelse(ifdouble,1,d_,) ppm98_interpolation_factors ^ (dxlag(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 c Get the advection fluxes c PPM98_MAP_FLUX_GAMMA_1D produces the advected zone interface fluxes in the c longitudinal direction (the direction of this pass) needed to c the new zone-averaged quantities on the grid defined by XXLLAG backc to the grid defined by XL, in 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 constructin c interpolations of the 0 Riemann invarient (adabatic) c (Returns UNSMAD) c 3) Using these piecewise parabolic interpolations for c , density, velocity and pressure, find the c averages in the appropriate overlap regions between c the two grids necessary to map the zone-averages c from XLLAG to XL. c Also, determine the remap courant number, which is c the largest fractional zone (in volume or mass) c advected. c (Returns DVOLL, DMASSL, DMOML, and DENL, the c fractions DVOLFL and DMASFL, and COURNO.) c c These individual steps are also available as separate c modules. c c map_flux gamma requires 4 boundary zones nbdy = 4 k=kstart2 nzones = kdo2-2*nbdy c c Note a 1d routine is correct, as the transverse velocity will be treated as c a passively advected quantity. C call ifelse(ifdouble,1,d_,,) ppm98_map_flux_gamma_1d ^ ( dalfac(K), darfac(K), facda(K), fcdazl(K), facdal(K), fa6dal(K), ^ fa6dar(K), flldal(K), flldar(K), frrdal(K), frrdar(K), ^ xllag(k), xl(K), rholag(k), plag(k), ulag(k), ^ unsmad(k), unsmup(k), stpadb(k), dvoll(k), ^ dmassl(k), dmoml(k), denl(k), dvolfl(k), dmassfl(k), 6 courno(k), dt, eosgam, smalle, nzones, nbdy) kstart3=kstart2+nbdy kend3=kend2-nbdy kdo3 = kdo2-2*nbdy 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 courmxa = courno(KSTART3 -1) do i=KSTART3 , KEND3+1 courmxa = max (courmxa, courno(i)) enddo c c These can be examined by c c print *,' i,rholag(i), dvoll(i), dmassl(i), dmoml(i), denl(i)' c do i=kstart3,kend3+1 c write(6,66) i,rholag(i), dvoll(i), dmassl(i), dmoml(i), denl(i) c enddo c c c New update the zone averages by conservatively differencing these c fluxes. Compute a new pressure if one wishes. c c print *,' i,rhonu(i), unu(i), eenew, pnu(i) ' do i=kstart3,kend3 dminv(i) = 1.0/dm(i) dmnu(i) = dm(i) + dmassl(i)- dmassl(i+1) dmnuinv(i) = 1.0/dmnu(i) rhonu(i) = dmnu(i) * dxinv(i) unu(i) = (ulag(i)*dm(i) + dmoml(i)-dmoml(i+1))/dmnu(i) eenew = (e_totlag(i)*dm(i) + denl(i)-denl(i+1))/dmnu(i) eknew = 0.5*unu(i)*unu(i) einew = max (smalle, eenew - eknew) c c NOTE this e_totnu does not yet contain a contribution from UT c e_totnu(i) = einew+eknew pnu(i) = (eosgam - 1.0) * rhonu(i)*einew c write(6,66) i,rhonu(i), unu(i), eenew, pnu(i) enddo c c Now, advect the transverse velocity ut: c c First, compute the advected volume fractions. Historicaly, PPM c has used a spatial fraction: dxl/dx, when logically it would c seem that a mass fraction should be used. I got no satisfactory c explaination from Woodward, as usual, other than "that's how I did c it." For now, I follow this tradition. I expect to change it. I c also agree with Woodward that the difference will probably be c negligable. I like consistancy. And in this formulation, advected c mass fractions are no more expensive than spatial or volume fractions. c For know, and in Cartesian coordinates (spatial = volume) we use c dvoll/dvol c do i=KSTART2,KEND2 c if( dvoll(i) .le. 0.0 ) then c dvolfl(i) = -dvoll(i) * dvolinv(i) c else c dvolfl(i) = dvoll(i) * dvolinv(i-1) c endif c dvolfl(i) = min (one, dvolfl(i)) c enddo c c Now passively advect the transverse volocity, using mass weightings: c c c c PPM98_ADVECT_PASSIVE0 performs the entire advection step. Its c c steps are: c 1) Interpolate, with discontinuity detection and c monotonicity constraints. (Returns UNSMUTUT c and STEEPUT) c 2) Integrate the piecewise parabolic interpolations c over the appropriate fraction of each zone c indicated by DVOLFL. (As noted above, logically c it would appear that this should be DMASFL, the c mass fraction.) c 3) Construct weighted fluxes from step 2) using the c provided interface weightings, here dmassl, c the mass advected at the lefthand zone interface. c (Returns DMOMTL) 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 volume (DMNUINV) as weightings. 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 c c Thus the following routine requires various kinds of input: c 1) zone averages of UT for interpolation. These will c require as many boundary cells as were required for c the call to PPM98_DEUL0_FLUX_GAMMA. c 2) Flux factors DVOLFL, DMASSL, required on a smaller c range than the zone averages in 1. This is because c they will be applied only at interfaces where valid c interpolations exist to be integrated over the indicated c advected zone fractions, and the interpolation routines c require 3 cells usually. The number of interface boundary c zone is passed here as a check. At this point, it should c usually be 1. c 3) the zone update weightings (DM, DMNUINV=1/DMNU), which c require no zone boundary zones. They need be valid only c on the range 1:nzones, the same as the updated c transverse velocity UTNU. c c All these arrays begin with the same offset K, for sanity. This c ensures that when refferencing array elements of "I", all such c references are to the same physical zone. Chaos wold result otherwise. 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 = KDO2 -2*nbdy_za K = KSTART2 call ifelse(ifdouble,1,d_,) ppm98_advect_passive0 ^ (dvolfl(K), dmassl(K), dm(K), ^ dmnuinv(K), ut(K), smalldu(K), ^ unsmutut(K), steeput(K), dmomtl(K), ^ utnu(K), nzones, nbdy_za ) 96 format(1h ,i4,1p,9(1x,e12.5)) c print *,' i,ut(i), steeput(i), unsmutut(i)' c do i=1,nzones c write(6,96) i,ut(i), steeput(i), unsmutut(i) c enddo ccc ccc Update the transverse kinetic energy and add it to the new total E do i=KSTART3 ,KEND3+1 if (dmassl(i) .eq. 0.0) then dektl(i) = 0.0 else dektl(i) = 0.5 * dmomtl(I)*dmomtl(I) / dmassl(i) endif enddo do i=KSTART3 ,KEND3 ekt = 0.5 * ut(I)*ut(I) ektnu = ( ekt * dm(i) + dektl(i) - dektl(i+1))*dmnuinv(i) e_totnu(i) = e_totnu(i) + ektnu enddo c C ********************************************************** C * * C * Done with the hydro step, begin diffusion step * C * * C ********************************************************** 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 KSTART3 C and the number of zone to process is KDO3 C c BUT REALIZE THAT PPM98_EUL_DIFFUSE_1D OPERATES ON THE INITIAL C DATA. c C KSTARTDF = KSTART3 KENDDF = KEND3 KDODF = KDO3 c c c Get difusion velocities c c For efficiency, only those zones in which the velocity interpolations c are "unsmooth" will be checked for shocks. Therefore, we combine the c results for the transvers velocity UT and longitudinal velcity u. We c are interested only in non zero, not actual values, se we can just add. do i=KSTART3 ,KEND3 unsmooth(i) = unsmutut(i) + unsmup(i) enddo c c PPM98_LAG_DIFUSE0_GAMMA_2D obtains zone centered Eulerian diffusion c velocities. c (It's like sausage, you don't really want to know it's done) c K=KSTARTDF nzones = KDODF c it is a uniform grid so any dx will do: dtbydx = dt * dxinv(5) call ifelse(ifdouble,1,d_,) ppm98_lag_difuse0_gamma_2d ^ (uzbb(k), u(k-2), u(k-1), u(K), u(k+1), u(k+2), ^ uztt(k), utzbb(k), utzb(k), ut(k-2), ut(K), ut(k+2), ^ utzt(k), utztt(k), pzbb(k), p(k-2), p(K), p(k+2), ^ pztt(k), rhozbb(k), rho(k-2), rho(K), rho(k+2), rhoztt(k), ^ ceul(K), unsmooth(k), difuse_vel(K), shockd(k), courno(K), ^ eosgam, dtbydx, smallu, nzones) c c update the zone counts KSTARTDF1 = KSTARTDF KENDDF1 = KENDDF KDODF1 = KDODF courdf1 = courno(KSTARTDF1) do i=KSTARTDF1+1,KENDDF1 courdf1 = max (courdf1,courno(i)) enddo c print *,' i,rho(i),difuse_vel(i),courno(i)' c do i=KSTARTDF1,KENDDF1 c write(6,66) i,rho(i),difuse_vel(i),courno(i) c enddo c At this point, we should have one extra difuseion velocity c at each end of the working arry. 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 First 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 ) c get the transverse momentum diffusion fluxes c c PPM98_DIFUSEFLUX_PASSIVE uses the the leftward and rightward c diffusion volumes (MASLL, MASRL) computed above to compute c lefthand zone interface mass weighted interface difusion momenta dmomtl. nbdy = 1 call ifelse(ifdouble,1,d_,) ppm98_difuseflux_passive ^ (utnu(K), dmasll(K), dmasrl(K), dmomtl(K), nzones, nbdy) 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 Now, 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 this saves a flop or two: if ((dvoll(i) + dvoll(i+1)) .ne. zero) then dmnew(i) = dmnu(i) + dmassl(i) - dmassl(i+1) rhonu(i) = dmnew(i) * dvolinv(i) dmnewi = one / dmnew(i) unu(i) = (unu(i) * dmnu(i) + dmoml(i) - dmoml(i+1)) * dmnewi utnu(i) = (utnu(i) * dmnu(i) + dmomtl(i) - dmomtl(i+1)) *dmnewi e_totnu(i) = (e_totnu(i) * dmnu(i) + ^ detotl(i) - detotl(i+1)) * dmnewi 9099 continue endif 9101 continue 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 utsq = utnu(i) * utnu(i) if ( utsq .lt. ssmlu2) then utnu(i) = 0.0 utsq = 0.0 endif eknu = 0.5 * (usq + utsq ) e_intnu = max (smalle, (e_totnu(i) - eknu)) pnu(i) = gamma1 * rhonu(i) * e_intnu 9100 continue c c The final Courant number is the maximum of the c hydrodynamics advective, and diffusive Courant numbers. c courmx = max (courmxh, courmxa, courdf1, courdf2 ) deallocate (smalldu, steeput, dmomtl, shockd) deallocate (dektl, dmnuinv, unsmutut, unsmooth) deallocate (courno, difuse_vel) deallocate ( unsmad, unsmup, stpadb, ceul) deallocate ( rhoavl, pavl, uavl) deallocate ( e_total, e_totnu, dx, dxinv) deallocate ( difusl, dmassl, dmoml, detotl) deallocate ( dmasll, dmasrl) deallocate ( xllag, dxlag) deallocate ( rholag, plag) deallocate ( e_totlag, ulag) deallocate (dvol, dvoll) deallocate (dm, dmnu, dmnew ) deallocate (dvolinv, dminv) deallocate (dmassfl, dvolfl, denl) c deallocate (dalfac, darfac, facda, fcdazl, ^ facdal, fa6dal, fa6dar, frrdal, ^ frrdar, flldal, flldar) c return end c