c c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** THIS CODE IS SPPM3D ***** c ***** THIS CODE WAS WRITTEN BY ***** c ***** ***** c ***** PAUL R. WOODWARD ***** c ***** ***** c ***** February, 1993 ***** c ***** REVISED March, 1994 ***** c ***** REVISED August, 1995 ***** c ***** Revised September, 1995 ***** c ***** ***** c ***** ***** c ***** COPYRIGHT Paul R. Woodward ***** c ***** March, 1994 ***** c ***** September, 1995 ***** c ***** ALL RIGHTS RESERVED ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c Laboratory for Computational Science & Engineering c University of Minnesota c c ================================================================ c iq.h c Defines per node tile size for the computational mesh. c c Following settings are an example of a 384-cubed total c problem size with a 2x2x1 nodelayout. c IQ should be set to at least the max of 2048, 2048, & 2048 c REAL Either "float" or "double", for desired fl. pt. size c For the proper result, also use a compiler flag to c change the default size of the 'real' declaration. c****************************************************** c switch settings to control alternate code in sppm.m4 c****************************************************** c c c======================================================================= c c c======================================================================= c c c======================================================================= c c THIS CODE PACKAGE IS SET UP ASSUMING THAT IT WILL FIRST BE c PREPROCESSED USING THE M4 UNIX MACRO PROCESSOR, BASED UPON c VALUES SET FOR THE ABOVE FLAGS AND CONSTANTS. c M4 MUST BE GIVEN THE OPTION -B30000 IN ORDER TO INCREASE c ITS BUFFER SIZES SUFFICIENTLY TO HANDLE THIS CODE. c c======================================================================= c c c c======================================================================= c c The flag 0 indicates whether or not the computation is to c be performed on a Cray computer. When this flag is 1, a Cray c computation is assumed, with 64-bit single precision arithmetic c and the special ISMAX intrinsic function for finding the maximum c value of a vector. When this flag is 0, a microprocessor-based c computation is assumed, in which double-precision 64-bit c arithmetic and associated intrinsic functions may be required. c The maximum Courant number will be determined from a simple loop c containing a test. c c======================================================================= c c The flag ifsngl indicates whether or not the computation is to c be performed using single precision. When this flag is 1, c single precision is used, when the flag is 0, double precision c is used. c ifsngl = 1 c c======================================================================= c c The flag ifcvmg indicates whether or not the intrinsic functions c cvmgm and cvmgz will be used to perform vectorizable logic. c On a Cray, this results in the fastest code. On a CM-5 Connection c machine, these intrinsic functions will be converted by the c Fortran-P precompiler to equivalent intrinsic functions, yielding c once again the fastest implementation. To use these special c vector logic functions, set the flag to 1. When this flag is set c to 0, equivalent IF constructions will be substituted. These c give the best performance on standard microprocessors. The Cray c compiler can also handle these constructions, but optimal code c will not result. c Also see the flag 1 described below. c ifcvmg = 0 c c======================================================================= c c The flag ifinln indicates whether or not the function cvmgm c will be "inlined." This flag takes effect only when the flag c 0 is set to 1. If 1 is set to 1, then all c occurences of cvmgm will be converted to an inlined version of c this function which will be recognized by the compiler as a c "conditional move" instruction. This will permit software c pipelining of loops containing cvmgm. That is, it will permit c microprocessor compilers to recognize that these loops containing c cvmgm can be vectorized. It may also provide increased c performance on microprocessors such as the DEC Alpha, which have c conditional move instructions implemented in hardware. c ifinln = 1 c c======================================================================= c c The flag 1 indicates whether or not square roots will be c taken in the normal way. If 1 is set to 1, then we will c compute the square root of a quantity a as sqrt (a). c If 1 is set to 0, then we will compute this square root c via a / sqrt (a) c This longer apparent computation can exploit special, fast vendor c implementations of the oft-needed reciprocal square root function. c On the Silicon Graphics MIPS R-8000 processor, the fast square c root computation can be obtained by setting 1 to 1 and c by invoking the compiler option -OPT:fast_sqrt c c======================================================================= c c The flag 1 indicates whether or not the computation will c be performed on a superscalar cache-based microprocessor with dual c sets of floating point units (like the IBM RIOS-2 and the MIPS c R-8000 chip sets). If such a computing platform will not be used, c 1 should be set to 0. When instead 1 is set to 1, c long vectorizable loops will be broken into shorter segments for c "software pipelining." This acknowledges the primitive state of c software pipelining technology in microprocessor compilers today. c Breaking these loops apart causes unnecessary memory traffic, but c on these microprocessors this memory traffic is restricted to the c cache, and therefore does not cause the performance degradation c which 1. might expect (although it is still bad, and these c microprocessor compilers should be improved). c Note this setting is also appropriate for vectorizing compilers. c ifsupr = 1 c c======================================================================= c c The constant IQ gives the maximum length of a physical row of c zones (not counting fake zones used to implement boundary condi- c tions) which will occur. For a given problem, it is generally c set to the maximum of the physical row lengths in the x-, y- and c z-directions. If arrays with larger rows are sent to this code c package by a calling routine, an error will result and the job c will be killed. c IQ = 2048 c c======================================================================= c c c c c c c c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** THIS CODE IS SPPM3D ***** c ***** THIS CODE WAS WRITTEN BY ***** c ***** ***** c ***** PAUL R. WOODWARD ***** c ***** ***** c ***** February, 1993 ***** c ***** REVISED March, 1994 ***** c ***** REVISED August, 1995 ***** c ***** Revised September, 1995 ***** c ***** ***** c ***** ***** c ***** COPYRIGHT Paul R. Woodward ***** c ***** March, 1994 ***** c ***** September, 1995 ***** c ***** ALL RIGHTS RESERVED ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c c c c c subroutine sppm (vvv1, vvv2, rhozbb, rhoztt, rhoznn, rhozff, 1 vvv3, pzbb, pztt, pznn, pzff, 2 vvv4, uxzbb, uxztt, uxznn, uxzff, 3 vvv5, uyzb, uyzt, uyzbb, uyztt, 4 uyznn, uyzff, 5 vvv6, uzzbb, uzztt, 6 uzzn, uzzf, uzznn, uzzff, 7 rhonu, pnu, uxnu, uynu, uznu, 8 dt, courmx, gamma, 9 smlrho, smalle, smallp, smallu, iii1, iii2) c c subroutine sppm (vvv1,vvv2,vvv3,vvv4,vvv5,vvv6,vvv7,vvv8,vvv9, c 1 vvv10,vvv11,vvv12,vvv13,vvv14,vvv15,vvv16,vvv17, c 2 vvv18,vvv19,vvv20,vvv21,vvv22,vvv23,vvv24,vvv25, c 3 vvv26,vvv27,vvv28,vvv29,vvv30,vvv31,vvv32,vvv33, c 4 vvv34,vvv35, c 5 ccc1,ccc2,ccc3,ccc4,ccc5,ccc6,ccc7, iii1,iii2) c c call sppm (xl, rho, rhozbb, rhoztt, rhoznn, rhozff, c 1 p, pzbb, pztt, pznn, pzff, c 2 ux, uxzbb, uxztt, uxznn, uxzff, c 3 uy, uyzb, uyzt, uyzbb, uyztt, c 4 uyznn, uyzff, c 5 uz, uzzbb, uzztt, c 6 uzzn, uzzf, uzznn, uzzff, c 7 rhonu, pnu, uxnu, uynu, uznu, c 8 dt, courmx, gamma, c 9 smlrho, smalle, smallp, smallu, nx, nbdy) c c We have a problem here with notation. We have consistently c used l and r for left and right, and also we have c used b and t for bottom and top. Now we need c designators for a third dimension. Up and down conflict c in meaning with bottom and top. Front and back are fine, but c they would introduce a second "b" designator. Therefore I c have chosen "near" and "far," with designators n and f, c to correspond to bottom and top in the other transverse c dimension. c c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** ***** c ***** ssss PPPPP ppppp m m ***** c ***** s s p p p p mm mm ***** c ***** s ppppp ppppp m m m m ***** c ***** s p p m m m ***** c ***** s s p p m m ***** c ***** ssss p p m m ***** c ***** ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c c parameter (nbdy=5) dimension vvv1(-nbdy+1:iii1+nbdy+1), vvv2(-nbdy+1:iii1+nbdy), 1 vvv3(-nbdy+1:iii1+nbdy), vvv4(-nbdy+1:iii1+nbdy), 2 vvv5(-nbdy+1:iii1+nbdy), vvv6(-nbdy+1:iii1+nbdy) c 3 vvv7(-nbdy+1:iii1+nbdy+1), vvv8(-nbdy+1:iii1+nbdy+1), c 4 vvv9(-nbdy+1:iii1+nbdy+1), vvv10(-nbdy+1:iii1+nbdy+1) c dimension vvv11(-nbdy+1:iii1+nbdy+1), vvv12(-nbdy+1:iii1+nbdy+1), c 1 vvv13(-nbdy+1:iii1+nbdy+1), vvv14(-nbdy+1:iii1+nbdy+1), c 2 vvv15(-nbdy+1:iii1+nbdy+1), vvv16(-nbdy+1:iii1+nbdy+1), c 3 vvv17(-nbdy+1:iii1+nbdy+1), vvv18(-nbdy+1:iii1+nbdy+1), c 4 vvv19(-nbdy+1:iii1+nbdy+1), vvv20(-nbdy+1:iii1+nbdy+1) c dimension vvv21(-nbdy+1:iii1+nbdy+1), vvv22(-nbdy+1:iii1+nbdy+1), c 1 vvv23(-nbdy+1:iii1+nbdy+1), vvv24(-nbdy+1:iii1+nbdy+1), c 2 vvv25(-nbdy+1:iii1+nbdy+1), vvv26(-nbdy+1:iii1+nbdy+1), c 3 vvv27(-nbdy+1:iii1+nbdy+1), vvv28(-nbdy+1:iii1+nbdy+1), c 4 vvv29(-nbdy+1:iii1+nbdy+1), vvv30(-nbdy+1:iii1+nbdy+1) c dimension vvv31(-nbdy+1:iii1+nbdy+1), vvv32(-nbdy+1:iii1+nbdy+1), c 1 vvv33(-nbdy+1:iii1+nbdy+1), vvv34(-nbdy+1:iii1+nbdy+1), c 2 vvv35(-nbdy+1:iii1+nbdy+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 March, 1993 ***** c ***** REVISED August, 1995 ***** c ***** Revised September, 1995 ***** c ***** ***** c ***** ***** c ***** COPYRIGHT Paul R. Woodward ***** c ***** March, 1994 ***** c ***** September, 1995 ***** c ***** ALL RIGHTS RESERVED ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c c c c c This routine performs the hydrodynamic timestep. c c call sppm (xl, rho, rhozbb, rhoztt, rhoznn, rhozff, c 1 p, pzbb, pztt, pznn, pzff, c 2 ux, uxzbb, uxztt, uxznn, uxzff, c 3 uy, uyzb, uyzt, uyzbb, uyztt, c 4 uyznn, uyzff, c 5 uz, uzzbb, uzztt, c 6 uzzn, uzzf, uzznn, uzzff, c 7 rhonu, pnu, uxnu, uynu, uznu, c 8 dt, courmx, gamma, c 9 smlrho, smalle, smallp, smallu, nx, nbdy) c c c c INPUT: c c c nbdy = the number of fake zones needed to implement boundary c conditions so that in the hydrodynamical calculation c all zones can be treated as if they were interior zones. c n = the number of real zones in the row to be processed. c dt = time step. c gamma = the ratio of specific heats. c smlrho = a trivial value of density. c smalle = a trivial value of energy per unit mass. c smallp = a trivial value of pressure. c smallu = a trivial value of velocity. c courmx = the maximum Courant number encountered so far (C*dt/dm, c or the largest fraction of a zone traversed by a c sound wave during the timestep). c c c The following arrays are given. Each allows space c for nbdy fake zones on either end of the physical row. Data c is not provided for these fake zones. c c xl = the locations of the left-hand interfaces of the zones. c rho = the average densities within the zones. c p = the average pressures within the zones. c ux = the zone-averaged longitudinal velocities. c uy = the zone-averaged transverse velocities. c uz = the zone-averaged transverse velocities in the third c dimension. c c Additional row arrays are provided which give values of these c physical variables in rows of the grid positioned 2 rows below c the row of interest in the y-direction (suffix "zbb") and 2 c rows above in the y-direction (suffix "ztt"), as well as 2 rows c nearer than the row of interest in the z-direction (suffix "znn") c and 2 rows further in the z-direction (suffix "zff"). c For the y- and z-velocities only, additional row arrays are c provided which give values of these physical variables in rows c of the grid positioned 1 row belowthe row of interest in the c y-direction (suffix "zb") and 1 row above in the y-direction c (suffix "zt"), as well as 1 row nearer than the row of interest c in the z-direction (suffix "zn") and 1 row further in the c z-direction (suffix "zf"). c c Note that the suffix "zb" denotes "zone on the bottom," c the suffix "zt" denotes "zone on the top," c the suffix "zn" denotes "zone nearer," c the suffix "zf" denotes "zone further," c and stuttering of the final letter in these suffixes denotes c 2 zones away rather than 1. c c The computation here will account only for derivatives in the c x-direction, the longitudinal direction for this "pass" through c the grid. Derivatives in the transverse directions, y and z, c will be used only in the routine which computes a diffusion c velocity. These derivatives will be used to sense the presence c and direction of propagation of shocks. c c c The model equation of state is: c c p = (gamma - 1) * rho * ei c c c OUTPUT: c c c courmx = the maximum Courant number encountered so far (C*dt/dm, c or the largest fraction of a zone traversed by a c sound wave during the timestep). c c c The following arrays are returned. Each allows space c for 2. fake zones on either end of the physical row. Data c need not be provided for these fake zones. However, xnul c must be set for the first fake zone. c c rhonu = new average densities within the zones. c pnu = new zone-averaged pressures. c uxnu = new zone-averaged longitudinal velocities. c uynu = new zone-averaged transverse velocities (y-components). c uznu = new zone-averaged transverse velocities (z-components). c c c c c c c c This routine performs SPPM hydrodynmaics in Lagrangian style c using a Riemann solver which treats all waves as discontinuities. c A simple gamma-law equation of state is used, and an initially c uniform grid with either periodic or continuation boundary c conditions is assumed. After the Lagrangian step, an SPPM remap c step is performed. c c SPPM is a simplified (hence the "S") version of PPM, the c Piecewise-Parabolic method. PPM was originally developed by c Paul R. Woodward and Philip Colella at the Lawrence Livermore c National Laboratory, extending work of Bram van Leer and Paul c Woodward on the earlier MUSCL hydrodynamics scheme. c SPPM contains several features of PPM, but significant features c have been omitted. Among these omitted features are contact c discontinuity detection and steepening and the computation of c a coefficient of numerical viscosity which adjusts to the local c needs of the flow in both space and time. SPPM is not intended c to be robust or general. Instead it is intended to be relatively c simple. Because it will be used for benchmarking purposes, it c contains a nonlinear Riemann solver and a careful computation of c the Courant time step limit. These computations serve to give c SPPM something more like the balance of computation and memory c traffic of the full PPM scheme. These computations are not c useless, although they certainly are not required for all c problems. c c parameter (iq=2048) parameter (jq=iq+nbdy) parameter (mjq=-nbdy+1) dimension xl(mjq:jq+1), rho(mjq:jq), p(mjq:jq), 1 ux(mjq:jq), uy(mjq:jq), uz(mjq:jq), 2 rhonu(mjq:jq), pnu(mjq:jq), uxnu(mjq:jq), 3 uynu(mjq:jq), uznu(mjq:jq), dxnui(mjq:jq) dimension rhozbb(mjq:jq), pzbb(mjq:jq), uxzbb(mjq:jq), 1 rhoztt(mjq:jq), pztt(mjq:jq), uxztt(mjq:jq), 2 rhoznn(mjq:jq), pznn(mjq:jq), uxznn(mjq:jq), 3 rhozff(mjq:jq), pzff(mjq:jq), uxzff(mjq:jq) dimension uyzbb(mjq:jq), uyzb(mjq:jq), uyzt(mjq:jq), 1 uyztt(mjq:jq), uyznn(mjq:jq), uyzff(mjq:jq), 2 uzzbb(mjq:jq), uzzn(mjq:jq), uzzf(mjq:jq), 3 uzztt(mjq:jq), uzznn(mjq:jq), uzzff(mjq:jq) dimension courno(mjq:jq), c(mjq:jq), v(mjq:jq), 1 uxpavl(mjq:jq), uxavl(mjq:jq), pavl(mjq:jq), 2 dux(mjq:jq), dp(mjq:jq), xf(mjq:jq), 3 xf1(mjq:jq), xnul(mjq:jq), dx(mjq:jq), 4 dm(mjq:jq), dtbydm(mjq:jq), vnu(mjq:jq), 5 dmnu(mjq:jq), dtbydx(mjq:jq), rho6(mjq:jq) dimension ux6(mjq:jq), p6(mjq:jq), uxl(mjq:jq), 1 uxr(mjq:jq), pl(mjq:jq), pr(mjq:jq), 2 dxnu(mjq:jq), drho(mjq:jq), dmassl(mjq:jq), 3 dmomxl(mjq:jq), dmomyl(mjq:jq), dmomzl(mjq:jq), 4 denl(mjq:jq), e(mjq:jq), enu(mjq:jq) dimension rhol(mjq:jq), rhor(mjq:jq), uy6(mjq:jq), 1 uyl(mjq:jq), uyr(mjq:jq), duy(mjq:jq), 2 uzl(mjq:jq), uzr(mjq:jq), duz(mjq:jq), 3 rplusl(mjq:jq), rplusr(mjq:jq), drplsl(mjq:jq), 4 drmnsl(mjq:jq), drplus(mjq:jq), drmnus(mjq:jq) dimension rmnusl(mjq:jq), rmnusr(mjq:jq), difuse(mjq:jq), 1 rplus6(mjq:jq), rmnus6(mjq:jq), uz6(mjq:jq), 2 dvoll(mjq:jq) c c dimension pll(mjq:jq), prl(mjq:jq), uuxavl(mjq:jq), 1 uxll(mjq:jq), uxrl(mjq:jq), dvolfl(mjq:jq), 2 dmfl(mjq:jq), dmfl1(mjq:jq), dvlfl1(mjq:jq) c c forthd = 4. / 3. small = 1.0e-08 c c 8 dt, courmx, gamma, c 9 smlrho, smalle, smallp, smallu, iii1, iii2) c write(0,*) "sppm: dt,courmx,gamma=", dt,courmx,gamma c write(0,*) "sppm: smlrho, smalle, smallp, smallu, iii1, iii2=", c & smlrho, smalle, smallp, smallu, iii1, iii2 if (iii2 .eq. nbdy) go to 10 write (6,11) iii2, nbdy 11 format (/'The value of NBDY passed to SPPM from the calling'/ 1 'routine is ',i3,' while that demanded by SPPM is ', 2 i3/'Job Killed.') stop 10 continue c if (iii1 .le. iq) go to 20 iiiq = iq write (6,21) iii1,iiiq 21 format (/'The value of N passed to SPPM from the calling'/ 1 'routine is ',i5,', while that allowed by the',/ 2 'dimension statements (parameter iq) in SPPM is only ', 3 i5/'Job Killed.') stop 20 continue c c dt = ccc1 c courmx = ccc2 c gamma = ccc3 c smlrho = ccc4 c smalle = ccc5 c smallp = ccc6 c smallu = ccc7 n = iii1 c c call sppm (xl, rho, rhozbb, rhoztt, rhoznn, rhozff, c 1 p, pzbb, pztt, pznn, pzff, c 2 ux, uxzbb, uxztt, uxznn, uxzff, c 3 uy, uyzb, uyzt, uyzbb, uyztt, c 4 uyznn, uyzff, c 5 uz, uzzbb, uzztt, c 6 uzzn, uzzf, uzznn, uzzff, c 7 rhonu, pnu, uxnu, uynu, uznu, c 8 dt, courmx, gamma, c 9 smlrho, smalle, smallp, smallu, nx, nbdy) c c call sppm (vvv1, vvv2, rhozbb, rhoztt, rhoznn, rhozff, c 1 vvv3, pzbb, pztt, pznn, pzff, c 2 vvv4, uxzbb, uxztt, uxznn, uxzff, c 3 vvv5, uyzb, uyzt, uyzbb, uyztt, c 4 uyznn, uyzff, c 5 vvv6, uzzbb, uzztt, c 6 uzzn, uzzf, uzznn, uzzff, c 7 rhonu, pnu, uxnu, uynu, uznu, c 8 dt, courmx, gamma, c 9 smlrho, smalle, smallp, smallu, iii1, iii2) c c On a Cray, the copying which is commented out in the loop below c is essentially a 0.-cost item. It guarantees freedom from c any effects of aliasing of the arguments of sppm. c However, on a microprocessor, with only a third of the memory c bandwidth per clock which is enjoyed by the Cray family of c machines, this loop costs more than the Riemann solver in the c computationally intensive section below in the Lagrangian step c of the hydrodynamics computation. Therefore, we will exploit c our knowledge that the adjacent strips are treated in a read-only c fashion and are used only in difuze in order to reduce the c cost of the algorithm as a whole by eliminating most of this c unnecessary copying. Actually, it can all be eliminated, but c the arrays xl, rho, p, and ux are rewritten during the c course of the calculation, and it is nice to have the arrays c rhonu, pnu, uxnu, uynu, and uznu kept separate from them c explicitly, in case the calling routine chooses to assign the c same storage to these quantities. c do 100 i = -nbdy+1,n+nbdy+1 xl(i) = vvv1(i) 100 continue do 110 i = -nbdy+1,n+nbdy rho(i) = vvv2(i) p(i) = vvv3(i) ux(i) = vvv4(i) uy(i) = vvv5(i) uz(i) = vvv6(i) 110 continue c rhozbb(i) = vvv3(i) c rhoztt(i) = vvv4(i) c rhoznn(i) = vvv5(i) c rhozff(i) = vvv6(i) c p(i) = vvv7(i) c pzbb(i) = vvv8(i) c pztt(i) = vvv9(i) c pznn(i) = vvv10(i) c pzff(i) = vvv11(i) c ux(i) = vvv12(i) c uxzbb(i) = vvv13(i) c uxztt(i) = vvv14(i) c uxznn(i) = vvv15(i) c uxzff(i) = vvv16(i) c uy(i) = vvv17(i) c uyzb(i) = vvv18(i) c uyzt(i) = vvv19(i) c uyzbb(i) = vvv20(i) c uyztt(i) = vvv21(i) c uyznn(i) = vvv22(i) c uyzff(i) = vvv23(i) c uz(i) = vvv24(i) c uzzbb(i) = vvv25(i) c uzztt(i) = vvv26(i) c uzzn(i) = vvv27(i) c uzzf(i) = vvv28(i) c uzznn(i) = vvv29(i) c uzzff(i) = vvv30(i) c 100 continue c c c c Set up essential quantities. c gaminv = 1. / gamma gamma1 = gamma - 1. gm1inv = 1. / gamma1 c do 1000 i = -nbdy+1,n+nbdy dx(i) = xl(i+1) - xl(i) dm(i) = rho(i) * dx(i) v(i) = 1. / rho(i) dtbydx(i) = dt / dx(i) dtbydm(i) = dtbydx(i) * v(i) c c(i) = sqrt (gamma * p(i) * rho(i)) c courno(i) = c(i) * dtbydm(i) xf(i) = .5 * min (1., courno(i)) xf1(i) = 1. - forthd * xf(i) 1000 continue c c c Now compute the diffusion speed, difuse. c This diffusion speed will be used to apply extra dissipation at c the end of the entire Eulerian hydrodynamics calculation. c noff = 4 call difuze (dtbydx, rho, rhozbb, rhoztt, rhoznn, rhozff, 1 p, pzbb, pztt, pznn, pzff, 2 ux, uxzbb, uxztt, uxznn, uxzff, 3 uy, uyzb, uyzt, uyzbb, uyztt, 4 uyznn, uyzff, 5 uz, uzzbb, uzztt, uzzn, uzzf, 6 uzznn, uzzff, 7 courno, difuse, gamma, smallu, n, nbdy, noff) c c c Compute the Riemann invariant differences, drplsl and c drmnsl, which are centered on the left-hand zone interfaces. c Use the equations: c dR+ = du + dp/C c dR- = du - dp/C c do 1100 i = -nbdy+2,n+nbdy thyng = 2. * (p(i) - p(i-1)) / (c(i) + c(i-1)) drplsl(i) = (ux(i) - ux(i-1)) + thyng drmnsl(i) = (ux(i) - ux(i-1)) - thyng 1100 continue c c c Obtain monotonized parabolae for the Riemann invariants. c Use the array pavl for scratch space. c noff = 1 call dintrf (rplusl,rplusr,drplus,rplus6,drplsl,pavl, 1 n,nbdy,noff) call dintrf (rmnusl,rmnusr,drmnus,rmnus6,drmnsl,pavl, 1 n,nbdy,noff) c c c Set up and solve Riemann problem for fluxes. c This set up involves obtaining the average values of the c Riemann invariants within the domains of dependence to the c left and to the right of the left-hand interface of the c zone. From these, the average pressures and velocities c in those domains can be estimated. These give the effective c left and right states for a Riemann problem, whose solution c gives an estimate of the time-averaged fluid state at the c left-hand interface during the time step. These time-averaged c states determine the fluxes used in the conservation laws to c update the zone-averaged specific volume and velocity (in a c succeeding loop). c Also compute new grid locations. c gammp1 = gamma + 1. hgamp1 = .5 * gammp1 c c c c csuperscalar code c do 3000 i = -nbdy+3,n+nbdy-1 c c FIRST COMPUTE AVERAGES OF VELOCITY AND PRESSURE IN THE c DOMAINS OF DEPENDENCE TO THE LEFT AND RIGHT OF THE LEFT-HAND c INTERFACE OF ZONE (I). IN CLASS, WE CALLED THESE DOMAIN c AVERAGES ULL AND PLL (FOR THE LEFT DOMAIN) AND URL AND PRL c (FOR THE RIGHT DOMAIN). c rpll = rplusr(i-1) - xf(i-1) * (drplus(i-1) 1 - xf1(i-1) * rplus6(i-1)) rmll = rmnusr(i-1) - xf(i-1) * (drmnus(i-1) 1 - xf1(i-1) * rmnus6(i-1)) uxll(i) = ux(i-1) + .5 * (rpll + rmll) pll(i) = p(i-1) + .5 * c(i-1) * (rpll - rmll) diffux = ux(i) - ux(i-1) c sux = sign (1., diffux) c breaks vectorization for stupid compilers sux = 1. if ( diffux .lt. 0 ) sux = -sux diffl = sux * (uxll(i) - ux(i-1)) diffr = sux * (ux(i) - uxll(i)) if (diffl .lt. 0.) uxll(i) = ux(i-1) if (diffr .lt. 0.) uxll(i) = ux(i) diffp = p(i) - p(i-1) c sp = sign (1., diffp) sp = 1. if ( diffp .lt. 0 ) sp= -sp diffl = sp * (pll(i) - p(i-1)) diffr = sp * (p(i) - pll(i)) if (diffl .lt. 0.) pll(i) = p(i-1) if (diffr .lt. 0.) pll(i) = p(i) rprl = rplusl(i) + xf(i) * (drplus(i) + xf1(i) * rplus6(i)) rmrl = rmnusl(i) + xf(i) * (drmnus(i) + xf1(i) * rmnus6(i)) uxrl(i) = ux(i) + .5 * (rprl + rmrl) prl(i) = p(i) + .5 * c(i) * (rprl - rmrl) diffl = sp * (prl(i) - p(i-1)) diffr = sp * (p(i) - prl(i)) if (diffl .lt. 0.) prl(i) = p(i-1) if (diffr .lt. 0.) prl(i) = p(i) diffl = sux * (uxrl(i) - ux(i-1)) diffr = sux * (ux(i) - uxrl(i)) if (diffl .lt. 0.) uxrl(i) = ux(i-1) if (diffr .lt. 0.) uxrl(i) = ux(i) 3000 continue c do 3100 i = -nbdy+3,n+nbdy-1 c c NOW, APPROXIMATING THE LAGRANGIAN SOUND SPEEDS IN THESE c DOMAINS OF DEPENDENCE BY THE AVERAGES OF THESE SPEEDS IN THE c ZONES TO THE LEFT AND RIGHT OF THE INTERFACE, COMPUTE A FIRST c GUESS AT THE TIME-AVERAGED PRESSURE, PAVL, AT THE INTERFACE c WHICH RESULTS FROM THE RIEMANN PROBLEM DESCRIBING THE INTERACTION c OF THE TWO DOMAINS. THIS FIRST GUESS IS OBTAINED BY ASSUMING c THAT THE TWO WAVES WHICH TRAVEL AWAY FROM THE INTERFACE ARE c WEAK. c dpavl = (pll(i) - prl(i) + c(i-1) * (uxll(i) - uxrl(i))) / 1 (c(i-1) + c(i)) pavl(i) = max (smallp, (prl(i) + dpavl * c(i))) uxavl(i) = uxrl(i) + dpavl c c NOW, FIND THE NONLINEAR WAVE SPEEDS WLL AND WRL FOR WAVES c WHICH RAISE (OR LOWER) THE PRESSURE TO THE VALUE PAVL JUST c ESTIMATED. USING THESE SPEEDS, NOW OBTAIN THE SLOPES OF THE c HUGONIOT CURVES FOR THE LEFT- AND RIGHT-HAND WAVES, DPDULL c AND DPDURL. ALSO OBTAIN THE VELOCITIES USTRLL AND USTRRL c WHICH LIE ON THESE HUGONIOT CURVES AT THE PRESSURE PAVL. c DPDULL AND DPDURL ARE THE SLOPES OF THE HUGONIOT CURVES AT c THE POINTS (USTRLL, PAVL) AND (USTRRL, PAVL). c wllfac = gamma1 * pll(i) + gammp1 * pavl(i) hrholl = .5 * rho(i-1) c wll = sqrt (hrholl * wllfac) c dpdull = wll * wllfac / (wllfac - hgamp1 * (pavl(i) - pll(i))) wrlfac = gamma1 * prl(i) + gammp1 * pavl(i) hrhorl = .5 * rho(i) c wrl = sqrt (hrhorl * wrlfac) c dpdurl = wrl * wrlfac / (wrlfac - hgamp1 * (pavl(i) - prl(i))) ustrll = uxll(i) - (pavl(i) - pll(i)) / wll ustrrl = uxrl(i) + (pavl(i) - prl(i)) / wrl c c NOW FIND THE POINT OF INTERSECTION OF THE TWO TANGENT LINES, c THE POINT (UAVL, PAVL), WHERE PAVL IS NOW AN IMPROVED c ESTIMATE OF THE TIME-AVERAGED PRESSURE AT THE INTERFACE AND c UAVL IS YOUR ESTIMATE OF THE TIME-AVERAGED VELOCITY THERE. c FOR SAFETY'S SAKE, YOU MAY WISH TO DEMAND THAT PAVL BE NO c LOWER THAN THE VALUE SMALLP, A TRIVIAL PRESSURE VALUE. c THIS WILL PREVENT THE GENERATION OF NEGATIVE PRESSURES IN c REGIONS OF FLOW CAVITATION. c thyng = (ustrrl - ustrll) * dpdurl / (dpdurl + dpdull) uxavl(i) = ustrll + thyng pavl(i) = max (smallp, (pavl(i) - thyng * dpdull)) 3100 continue c do 3200 i = -nbdy+3,n+nbdy-1 c c NOW, FIND THE NONLINEAR WAVE SPEEDS WLL AND WRL FOR WAVES c WHICH RAISE (OR LOWER) THE PRESSURE TO THE VALUE PAVL JUST c ESTIMATED. USING THESE SPEEDS, NOW OBTAIN THE SLOPES OF THE c HUGONIOT CURVES FOR THE LEFT- AND RIGHT-HAND WAVES, DPDULL c AND DPDURL. ALSO OBTAIN THE VELOCITIES USTRLL AND USTRRL c WHICH LIE ON THESE HUGONIOT CURVES AT THE PRESSURE PAVL. c DPDULL AND DPDURL ARE THE SLOPES OF THE HUGONIOT CURVES AT c THE POINTS (USTRLL, PAVL) AND (USTRRL, PAVL). c wllfac = gamma1 * pll(i) + gammp1 * pavl(i) hrholl = .5 * rho(i-1) c wll = sqrt (hrholl * wllfac) c dpdull = wll * wllfac / (wllfac - hgamp1 * (pavl(i) - pll(i))) wrlfac = gamma1 * prl(i) + gammp1 * pavl(i) hrhorl = .5 * rho(i) c wrl = sqrt (hrhorl * wrlfac) c dpdurl = wrl * wrlfac / (wrlfac - hgamp1 * (pavl(i) - prl(i))) ustrll = uxll(i) - (pavl(i) - pll(i)) / wll ustrrl = uxrl(i) + (pavl(i) - prl(i)) / wrl c c NOW FIND THE POINT OF INTERSECTION OF THE TWO TANGENT LINES, c THE POINT (UAVL, PAVL), WHERE PAVL IS NOW AN IMPROVED c ESTIMATE OF THE TIME-AVERAGED PRESSURE AT THE INTERFACE AND c UAVL IS YOUR ESTIMATE OF THE TIME-AVERAGED VELOCITY THERE. c FOR SAFETY'S SAKE, YOU MAY WISH TO DEMAND THAT PAVL BE NO c LOWER THAN THE VALUE SMALLP, A TRIVIAL PRESSURE VALUE. c THIS WILL PREVENT THE GENERATION OF NEGATIVE PRESSURES IN c REGIONS OF FLOW CAVITATION. c thyng = (ustrrl - ustrll) * dpdurl / (dpdurl + dpdull) uxavl(i) = ustrll + thyng pavl(i) = max (smallp, (pavl(i) - thyng * dpdull)) c c c NOW COMPUTE THE TOTAL ENERGY FLUX FROM THE PRODUCT OF UAVL c AND PAVL. YOU SHOULD ALSO COMPUTE THE NEW LOCATION OF THE c LEFT-HAND ZONE INTERFACE HERE. c uxpavl(i) = uxavl(i) * pavl(i) dvoll(i) = dt * uxavl(i) xnul(i) = xl(i) + dvoll(i) 3200 continue c csuperscalar code c c c c Apply conservation laws to obtain new zone-averaged c values of density and velocity. Also revise the c estimate of the Courant number to prevent zones from c tangling. c c Do not add the transverse kinetic energy term c .5 * (uy(i) * uy(i) + uz(i) * uz(i)) c to either enu(i) or to the expression in einu below. c This term just cancels out and requires additional computation. c c To make 32-bit arithmetic effective, we use formulae which c involve the smallest rounding errors, even if this costs us c 1. or 2. extra flops. c c do 4000 i = -nbdy+3,n+nbdy-2 c dxnu(i) = xnul(i+1) - xnul(i) dxnu(i) = dx(i) - dt * (uxavl(i) - uxavl(i+1)) courno(i) = max (courno(i), 1 (dtbydx(i) * (uxavl(i) - uxavl(i+1)))) c rhonu(i) = dm(i) / dxnu(i) vnu(i) = v(i) - dtbydm(i) * (uxavl(i) - uxavl(i+1)) rhonu(i) = 1. / vnu(i) c c This form of the conservation law involves more computation, c but it gives a better result for the new density when only c 32-bit arithmetic and inexact reciprocals are employed. c The reason for this is that the velocity difference will c vanish when it is supposed to, even for 32-bit arithmetic. c uxnu(i) = ux(i) + dtbydm(i) * (pavl(i) - pavl(i+1)) eee = gm1inv * p(i) * v(i) + .5 * ux(i) * ux(i) enu(i) = eee + dtbydm(i) * (uxpavl(i) - uxpavl(i+1)) einu = max (smalle, 1 (enu(i) - .5 * uxnu(i) * uxnu(i))) pnu(i) = gamma1 * rhonu(i) * einu 4000 continue c c c Now rename variables in preparation for the remap step. c No renaming is required for the transverse velocities uy c and uz. c do 5000 i = -nbdy+3,n+nbdy-2 rho(i) = rhonu(i) v(i) = vnu(i) p(i) = pnu(i) ux(i) = uxnu(i) dx(i) = dxnu(i) e(i) = enu(i) + 1 .5 * (uy(i) * uy(i) + uz(i) * uz(i)) 5000 continue c do 5100 i = -nbdy+3,n+nbdy-1 thyng = xl(i) xl(i) = xnul(i) xnul(i) = thyng 5100 continue c c c Obtain monotonized parabolae for rho, p, ux, uy, and uz. c Use the arrays uxavl and pavl for scratch space. c Ignore the fact that the grid is now nonuniform. c noff = 3 call interf (rho,rhol,rhor,drho,rho6,uxavl,pavl,n,nbdy,noff) call interf (p, pl, pr, dp, p6, uxavl, pavl, n, nbdy, noff) call interf (ux, uxl, uxr, dux, ux6, uxavl, pavl, n, nbdy, noff) call interf (uy, uyl, uyr, duy, uy6, uxavl, pavl, n, nbdy, noff) call interf (uz, uzl, uzr, duz, uz6, uxavl, pavl, n, nbdy, noff) c c c FROM THE NEW AND OLD X-COORDINATES OF THE LEFT-HAND ZONE c INTERFACES, COMPUTE DVOLL, THE VOLUME OF THE LAGRANGIAN c ZONE WHICH IS ADVECTED ACROSS THE LEFT-HAND INTERFACE OF c THE NEW EULERIAN ZONE. ADVECTION TO THE RIGHT GIVES A c POSITIVE ADVECTED VOLUME, WHILE ADVECTION TO THE LEFT RESULTS c IN A NEGATIVE ADVECTED VOLUME. c FROM THE DENSITY DISTRIBUTIONS DETERMINED ABOVE, COMPUTE THE c AVERAGE VALUES OF DENSITY, RHOAVL, IN THE VOLUMES ADVECTED c ACROSS THE LEFT-HAND ZONE INTERFACES. c FROM THESE AVERAGE DENSITIES AND THE ADVECTED VOLUMES, GET c THE ADVECTED MASSES, DMASSL. NOW USE THE INTERPOLATED c STRUCTURES FOR THE VELOCITIES, U AND UT, TO COMPUTE THE c AVERAGES OF THE VELOCITIES WITHIN THESE ADVECTED MASSES, AND c HENCE THE ADVECTED MOMENTA, DMOML AND DMOMTL. FROM THE c AVERAGES OF DENSITY, PRESSURE, AND THE VELOCITIES IN THE c ADVECTED MATERIAL, COMPUTE THE AVERGAES OF TOTAL ENERGY, E, c IN THE ADVECTED MASSES, AND HENCE THE ADVECTED TOTAL ENERGIES, c DENL. c c BE SURE TO RESET COURNO IF APPROPRIATE. c c FOR THE PROPER ADVECTION OF CONTACT DISCONTINUITIES, IT IS c ESSENTIAL THAT THE ADVECTED ENERGY, DENL, BE CONSTRUCTED c EXPLICITLY AS THE ADVECTED MASS TIMES HALF THE SQUARE OF THE c AVERAGE ADVECTED VELOCITY PLUS THE ADVECTED VOLUME TIMES c THE AVERAGE PRESSURE IN THE ADVECTED VOLUME TIMES THE c CONSTANT (1/(GAMMA-1)). THIS CONSTRUCTION GUARANTEES THAT c THE SIMPLE ADVECTION OF A CONTACT DISCONTINUITY GENERATES c NO SOUND WAVE DISTURBANCES (TO COMPUTER ROUND-OFF). c c twelth = .25 / 3. c c c c c csuperscalar code. Again modified for stupd vectorizing compiler c do 7000 i = -nbdy+5,n+nbdy-3 c dvoll(i) = (xl(i) - xnul(i)) thyng = dx(i-1) if (dvoll(i) .lt. 0.) then thyng = - dx(i) c else c thyng = dx(i-1) endif dvolfl(i) = dvoll(i) / thyng dvlfl1(i) = 1. - dvolfl(i) rrho = rho(i-1) rrhol = rhor(i-1) ddrho = -drho(i-1) if (dvoll(i) .lt. 0.) then rrho = rho(i) rrhol = rhol(i) ddrho = drho(i) c else c rrho = rho(i-1) c rrhol = rhor(i-1) c ddrho = - drho(i-1) endif thyng = rrho - rrhol rhoavl = rrhol + dvolfl(i) * (thyng 1 + dvlfl1(i) * (2. * thyng - ddrho)) dmassl(i) = dvoll(i) * rhoavl thyng = v(i-1) if (dvoll(i) .lt. 0.) then thyng = v(i) c else c thyng = v(i-1) endif dmfl(i) = dvolfl(i) * rhoavl * thyng dmfl1(i) = 1. - dmfl(i) courno(i) = max (courno(i), dvolfl(i), dmfl(i)) uux = ux(i-1) uuxl = uxr(i-1) ddux = - dux(i-1) if (dvoll(i) .lt. 0.) then uux = ux(i) uuxl = uxl(i) ddux = dux(i) c else c uux = ux(i-1) c uuxl = uxr(i-1) c ddux = - dux(i-1) endif thyng = uux - uuxl uuxavl(i) = uuxl + dmfl(i) * (thyng 1 + dmfl1(i) * (2. * thyng - ddux)) dmomxl(i) = dmassl(i) * uuxavl(i) 7000 continue c call dummy c do 7100 i = -nbdy+5,n+nbdy-3 uuy = uy(i-1) uuyl = uyr(i-1) dduy = - duy(i-1) if (dvoll(i) .lt. 0.) then uuy = uy(i) uuyl = uyl(i) dduy = duy(i) c else c uuy = uy(i-1) c uuyl = uyr(i-1) c dduy = - duy(i-1) endif thyng = uuy - uuyl uuyavl = uuyl + dmfl(i) * (thyng 1 + dmfl1(i) * (2. * thyng - dduy)) dmomyl(i) = dmassl(i) * uuyavl uuz = uz(i-1) uuzl = uzr(i-1) dduz = - duz(i-1) if (dvoll(i) .lt. 0.) then uuz = uz(i) uuzl = uzl(i) dduz = duz(i) c else c uuz = uz(i-1) c uuzl = uzr(i-1) c dduz = - duz(i-1) endif thyng = uuz - uuzl uuzavl = uuzl + dmfl(i) * (thyng 1 + dmfl1(i) * (2. * thyng - dduz)) dmomzl(i) = dmassl(i) * uuzavl pp = p(i-1) ppl = pr(i-1) ddp = - dp(i-1) if (dvoll(i) .lt. 0.) then pp = p(i) ppl = pl(i) ddp = dp(i) c else c pp = p(i-1) c ppl = pr(i-1) c ddp = - dp(i-1) endif thyng = pp - ppl ppavl = ppl + dvolfl(i) * (thyng 1 + dvlfl1(i) * (2. * thyng - ddp)) denl(i) = .5 * (uuxavl(i) * dmomxl(i) 1 + uuyavl * dmomyl(i) + uuzavl * dmomzl(i)) 2 + gm1inv * ppavl * dvoll(i) 7100 continue c csuperscalar code c c c c Now apply the conservation laws. c do 8000 i = -nbdy+5,n+nbdy-4 dxnu(i) = xnul(i+1) - xnul(i) dmnu(i) = dm(i) + dmassl(i) - dmassl(i+1) dxnui(i) = 1. / dxnu(i) rhonu(i) = dmnu(i) * dxnui(i) dmninv = 1. / dmnu(i) uxnu(i) = (ux(i) * dm(i) + dmomxl(i) - dmomxl(i+1)) * dmninv uynu(i) = (uy(i) * dm(i) + dmomyl(i) - dmomyl(i+1)) * dmninv uznu(i) = (uz(i) * dm(i) + dmomzl(i) - dmomzl(i+1)) * dmninv enu(i) = (e(i) * dm(i) + denl(i) - denl(i+1)) * dmninv 8000 continue c c c Now perform the diffusion step. This diffusion will matter c only in shocks which move slowly relative to the grid. c Without this diffusion step, these shocks would emit noise in c the Riemann invariant fields which cross the shock. c of impenetrable obstacle. c do 9000 i = -nbdy+6,n+nbdy-4 difusl = max (difuse(i-1), difuse(i)) dwoll = dt * difusl courno(i) = max (courno(i), ((dwoll + dwoll) * 1 max (dxnui(i-1), dxnui(i)))) ddmll = dwoll * rhonu(i-1) ddmrl = dwoll * rhonu(i) dmassl(i) = ddmll - ddmrl dmomxl(i) = ddmll * uxnu(i-1) - ddmrl * uxnu(i) dmomyl(i) = ddmll * uynu(i-1) - ddmrl * uynu(i) dmomzl(i) = ddmll * uznu(i-1) - ddmrl * uznu(i) denl(i) = ddmll * enu(i-1) - ddmrl * enu(i) 9000 continue c c do 9100 i = -nbdy+6,n+nbdy-5 dmnew = dmnu(i) + dmassl(i) - dmassl(i+1) rhonu(i) = dmnew * dxnui(i) dmnewi = 1. / dmnew uxnu(i) = (uxnu(i) * dmnu(i) + dmomxl(i) - dmomxl(i+1)) * dmnewi uynu(i) = (uynu(i) * dmnu(i) + dmomyl(i) - dmomyl(i+1)) * dmnewi uznu(i) = (uznu(i) * dmnu(i) + dmomzl(i) - dmomzl(i+1)) * dmnewi eknu = .5 * (uxnu(i) * uxnu(i) + uynu(i) * uynu(i) 1 + uznu(i) * uznu(i)) enew = (enu(i) * dmnu(i) + denl(i) - denl(i+1)) * dmnewi einu = max (smalle, (enew - eknu)) pnu(i) = gamma1 * rhonu(i) * einu 9100 continue c c c c Reset maximum Courant number, if appropriate. c c csuncode( c do 9200 i = 1,n courmx = max (courmx, courno(i)) 9200 continue c csuncode) c c call sppm (xl, rho, rhozbb, rhoztt, rhoznn, rhozff, c 1 p, pzbb, pztt, pznn, pzff, c 2 ux, uxzbb, uxztt, uxznn, uxzff, c 3 uy, uyzb, uyzt, uyzbb, uyztt, c 4 uyznn, uyzff, c 5 uz, uzzbb, uzztt, c 6 uzzn, uzzf, uzznn, uzzff, c 7 rhonu, pnu, uxnu, uynu, uznu, c 8 dt, courmx, gamma, c 9 smlrho, smalle, smallp, smallu, nx, nbdy) c c call sppm (vvv1, vvv2, rhozbb, rhoztt, rhoznn, rhozff, c 1 vvv3, pzbb, pztt, pznn, pzff, c 2 vvv4, uxzbb, uxztt, uxznn, uxzff, c 3 vvv5, uyzb, uyzt, uyzbb, uyztt, c 4 uyznn, uyzff, c 5 vvv6, uzzbb, uzztt, c 6 uzzn, uzzf, uzznn, uzzff, c 7 rhonu, pnu, uxnu, uynu, uznu, c 8 dt, courmx, gamma, c 9 smlrho, smalle, smallp, smallu, iii1, iii2) c c On a Cray, the copying which is commented out in the loop below c is essentially a 0.-cost item. It guarantees freedom from c any effects of aliasing of the arguments of sppm. c However, on a microprocessor, with only a third of the memory c bandwidth per clock which is enjoyed by the Cray family of c machines, this loop costs more than the Riemann solver in the c computationally intensive section below in the Lagrangian step c of the hydrodynamics computation. Therefore, we will exploit c our knowledge that the adjacent strips are treated in a read-only c fashion and are used only in difuze in order to reduce the c cost of the algorithm as a whole by eliminating most of this c unnecessary copying. Actually, it can all be eliminated, but c the arrays xl, rho, p, and ux are rewritten during the c course of the calculation, and it is nice to have the arrays c rhonu, pnu, uxnu, uynu, and uznu kept separate from them c explicitly, in case the calling routine chooses to assign the c same storage to these quantities. c c c ccc2 = courmx c c do 9900 i = -nbdy+6,n+nbdy-5 c vvv31(i) = rhonu(i) c vvv32(i) = pnu(i) c vvv33(i) = uxnu(i) c vvv34(i) = uynu(i) c vvv35(i) = uznu(i) c9900 continue c c return end c c---------------------------------------------------------------------72 c subroutine difuze (dtbydx, rho, rhozbb, rhoztt, rhoznn, rhozff, 1 p, pzbb, pztt, pznn, pzff, 2 ux, uxzbb, uxztt, uxznn, uxzff, 3 uy, uyzb, uyzt, uyzbb, uyztt, 4 uyznn, uyzff, 5 uz, uzzbb, uzztt, uzzn, uzzf, 6 uzznn, uzzff, 7 courno, difuse, gamma, smallu, n, nbdy, noff) c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** ***** c ***** dddd iii fffff u u zzzzz eeeee ***** c ***** d d i f u u z e ***** c ***** d d i fff u u z eee ***** c ***** d d i f u u z e ***** c ***** d d i f u u z e ***** c ***** dddd iii f uuu zzzzzz eeeee ***** c ***** ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c dimension ux(-nbdy+1:n+nbdy+1), uy(-nbdy+1:n+nbdy+1), & rho(-nbdy+1:n+nbdy+1), p(-nbdy+1:n+nbdy+1), & courno(-nbdy+1:n+nbdy+1), dtbydx(-nbdy+1:n+nbdy+1), & uxzbb(-nbdy+1:n+nbdy+1), uxztt(-nbdy+1:n+nbdy+1), & uyzbb(-nbdy+1:n+nbdy+1), uyztt(-nbdy+1:n+nbdy+1), & uyzb(-nbdy+1:n+nbdy+1), uyzt(-nbdy+1:n+nbdy+1), & uzzbb(-nbdy+1:n+nbdy+1), uzztt(-nbdy+1:n+nbdy+1), & pzbb(-nbdy+1:n+nbdy+1), pztt(-nbdy+1:n+nbdy+1), & rhozbb(-nbdy+1:n+nbdy+1), rhoztt(-nbdy+1:n+nbdy+1), & difuse(-nbdy+1:n+nbdy+1), uz(-nbdy+1:n+nbdy+1) dimension uxznn(-nbdy+1:n+nbdy+1), uxzff(-nbdy+1:n+nbdy+1), & uyznn(-nbdy+1:n+nbdy+1), uyzff(-nbdy+1:n+nbdy+1), & uzznn(-nbdy+1:n+nbdy+1), uzzff(-nbdy+1:n+nbdy+1), & uzzn(-nbdy+1:n+nbdy+1), uzzf(-nbdy+1:n+nbdy+1), & pznn(-nbdy+1:n+nbdy+1), pzff(-nbdy+1:n+nbdy+1), & rhoznn(-nbdy+1:n+nbdy+1), rhozff(-nbdy+1:n+nbdy+1) c common / diffuz / shkjmp, difcon, wavhaf c c parameter (nbdry=5) parameter (iq=2048) parameter (jq=iq+nbdry+1) parameter (mjq=-nbdry+1) dimension rho0x(mjq:jq), rho0y(mjq:jq), rho0z(mjq:jq), & rho1x(mjq:jq), rho1y(mjq:jq), rho1z(mjq:jq), & p0x(mjq:jq), p0y(mjq:jq), p0z(mjq:jq), & p1x(mjq:jq), p1y(mjq:jq), p1z(mjq:jq), & ux0x(mjq:jq), ux0y(mjq:jq), ux0z(mjq:jq), & ux1x(mjq:jq), ux1y(mjq:jq), ux1z(mjq:jq), & uy0x(mjq:jq), uy0y(mjq:jq), uy0z(mjq:jq), & uy1x(mjq:jq), uy1y(mjq:jq), uy1z(mjq:jq), & uz0x(mjq:jq), uz0y(mjq:jq), uz0z(mjq:jq), & uz1x(mjq:jq), uz1y(mjq:jq), uz1z(mjq:jq) dimension rhopre(mjq:jq), ppre(mjq:jq), dffuse(mjq:jq), & rhopst(mjq:jq), ppst(mjq:jq), & uxpre(mjq:jq), uypre(mjq:jq), uzpre(mjq:jq), & uxpst(mjq:jq), uypst(mjq:jq), uzpst(mjq:jq), & factrx(mjq:jq), factry(mjq:jq), factrz(mjq:jq) 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 March, 1993 ***** c ***** REVISED August, 1995 ***** c ***** Revised September, 1995 ***** c ***** ***** c ***** ***** c ***** COPYRIGHT Paul R. Woodward ***** c ***** March, 1994 ***** c ***** September, 1995 ***** c ***** ALL RIGHTS RESERVED ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c c c c This routine locates zones which contain the internal structure c of shocks and computes dissipation factors for these zones. c c c c Input to this routine: c c c n = the number of real, physical zones in the grid. c nbdy = the number of fake zones at each end of the physical c grid. c noff = the offset from each end of the full (real + fake) grid c which defines the strip in which meaningful values for c the ouput vector difuse is to be provided by this c routine. c c It is assumed that boundary conditions have been applied in the c calling program, so that all n+2*(nbdy-noff+2) zones of average c values are set properly. c c dt = the timestep. c gamma = the EOS is: p = (gamma-1)*rho*e_internal c smallu = a trivial positive value of velocity. c c c Row Arrays: c c dtbydx = the timestep divided by the zone widths. c rho = the zone-averaged densities. c rhozbb = the zone-averaged densities 2. rows below. c rhoztt = the zone-averaged densities 2. rows above. c rhoznn = the zone-averaged densities 2. rows nearer. c rhozff = the zone-averaged densities 2. rows further away. c p = the zone-averaged pressures. c pzbb = the zone-averaged pressures 2. rows below. c pztt = the zone-averaged pressures 2. rows above. c pznn = the zone-averaged pressures 2. rows nearer. c pzff = the zone-averaged pressures 2. rows further away. c ux = the zone-averaged x-velocities. c uxzbb = the zone-averaged x-velocities 2. rows below. c uxztt = the zone-averaged x-velocities 2. rows above. c uxznn = the zone-averaged x-velocities 2. rows nearer. c uxzff = the zone-averaged x-velocities 2. rows further away. c uy = the zone-averaged y-velocities. c uyzb = the zone-averaged y-velocities 1. row below. c uyzt = the zone-averaged y-velocities 1. row above. c uyzbb = the zone-averaged y-velocities 2. rows below. c uyztt = the zone-averaged y-velocities 2. rows above. c uyznn = the zone-averaged y-velocities 2. rows nearer. c uyzff = the zone-averaged y-velocities 2. rows further away. c uz = the zone-averaged z-velocities. c uzzbb = the zone-averaged z-velocities 2. rows below. c uzztt = the zone-averaged z-velocities 2. rows above. c uzzn = the zone-averaged z-velocities 1. row nearer. c uzzf = the zone-averaged z-velocities 1. row further away. c uzznn = the zone-averaged z-velocities 2. rows nearer. c uzzff = the zone-averaged z-velocities 2. rows further away. c courno = the Courant numbers in the zones. c c c c Output from this routine: c c c courno = Courant numbers, revised if necessary to account for c the Eulerian shock speeds computed here. c c difuse = the Eulerian diffusion velocity. c c difuse is a diffusion velocity which is zone-centered. It is c used in the remap step to prevent oscillations behind strong c shocks which are nearly aligned with the grid. To construct a c diffusion velocity at a zone interface in the remap step, we will c take the maximum of difuse in the zones on either side. c c c Dimensionless constants: c c c difcon = an overall factor multiplying the diffusion velocity c difuse. In the limit of a nearly steady strong shock, c difuse is difcon times the jump in the velocity. c It is recommended that difcon be of order unity. c A good value for difcon is 1/2. c shkjmp = the value of the fractional pressure jump, pjumpf, c below which we no longer consider a zone to be in or c near a shock. If a zone is not "shocked" we set the c diffusion velocity difuse to 0.. c The recommended value for shkjmp is 0.25, although c a value of 0.33 is also fine. c c small = 1.0e-08 difcon = .05 shkjmp = .25 c c c c c **************************************************** c **************************************************** c THE FORMULA BELOW FOR DIVU ASSUMES A UNIFORM GRID. c **************************************************** c **************************************************** c c difuse is used to reduce the generation of waves in Riemann c invariant fields which cross the shock. c c smlusq = smallu * smallu mbdy = nbdy - noff c c c c csuperscalar code c c do 3100 i = -mbdy+1,n+mbdy c c Now pick the pre- and post-shock values (0 and 1) along each c direction. c ddpx = p(i+2) - p(i-2) if (ddpx .lt. 0.) then p0x(i) = p(i+2) p1x(i) = p(i-2) ux0x(i) = ux(i+2) ux1x(i) = ux(i-2) uy0x(i) = uy(i+2) uy1x(i) = uy(i-2) uz0x(i) = uz(i+2) uz1x(i) = uz(i-2) rho0x(i) = rho(i+2) rho1x(i) = rho(i-2) else p0x(i) = p(i-2) p1x(i) = p(i+2) ux0x(i) = ux(i-2) ux1x(i) = ux(i+2) uy0x(i) = uy(i-2) uy1x(i) = uy(i+2) uz0x(i) = uz(i-2) uz1x(i) = uz(i+2) rho0x(i) = rho(i-2) rho1x(i) = rho(i+2) endif c 3100 continue c c c do 3200 i = -mbdy+1,n+mbdy c ddpy = pztt(i) - pzbb(i) if (ddpy .lt. 0.) then p0y(i) = pztt(i) p1y(i) = pzbb(i) ux0y(i) = uxztt(i) ux1y(i) = uxzbb(i) uy0y(i) = uyztt(i) uy1y(i) = uyzbb(i) uz0y(i) = uzztt(i) uz1y(i) = uzzbb(i) rho0y(i) = rhoztt(i) rho1y(i) = rhozbb(i) else p0y(i) = pzbb(i) p1y(i) = pztt(i) ux0y(i) = uxzbb(i) ux1y(i) = uxztt(i) uy0y(i) = uyzbb(i) uy1y(i) = uyztt(i) uz0y(i) = uzzbb(i) uz1y(i) = uzztt(i) rho0y(i) = rhozbb(i) rho1y(i) = rhoztt(i) endif c 3200 continue c c c do 3300 i = -mbdy+1,n+mbdy c ddpz = pzff(i) - pznn(i) if (ddpz .lt. 0.) then p0z(i) = pzff(i) p1z(i) = pznn(i) ux0z(i) = uxzff(i) ux1z(i) = uxznn(i) uy0z(i) = uyzff(i) uy1z(i) = uyznn(i) uz0z(i) = uzzff(i) uz1z(i) = uzznn(i) rho0z(i) = rhozff(i) rho1z(i) = rhoznn(i) else p0z(i) = pznn(i) p1z(i) = pzff(i) ux0z(i) = uxznn(i) ux1z(i) = uxzff(i) uy0z(i) = uyznn(i) uy1z(i) = uyzff(i) uz0z(i) = uzznn(i) uz1z(i) = uzzff(i) rho0z(i) = rhoznn(i) rho1z(i) = rhozff(i) endif 3300 continue c c ***************************************************************** c OPPS. I THINK THAT DEFINING "NEAR" TO BE LOWER VALUES OF Z c ***************************************************************** c PRODUCES A LEFT-HANDED COORDINATE SYSTEM. c ***************************************************************** c c do 3400 i = -mbdy+1,n+mbdy ddxux = ux(i+2) - ux(i-2) ddyuy = uyztt(i) - uyzbb(i) ddzuz = uzzff(i) - uzznn(i) dffuse(i) = difcon * abs (ddxux + ddyuy + ddzuz) ddxuxm = min (ddxux, 0.) ddyuym = min (ddyuy, 0.) ddzuzm = min (ddzuz, 0.) factrx(i) = ddxuxm * ddxuxm factry(i) = ddyuym * ddyuym factrz(i) = ddzuzm * ddzuzm deenom = 1. / (factrx(i) + factry(i) + 1 factrz(i) + smlusq) factry(i) = factry(i) * deenom factrz(i) = factrz(i) * deenom factrx(i) = 1. - (factry(i) + factrz(i)) 3400 continue c c c We will use the factors factor and factrt to blend data c from the longitudinal and transverse directions in order to c estimate shock jumps, orientation, speeds, etc., below. c This data blending technique generates more algebra here, c but it eliminates any tendency to flip-flop near a 45 degree c shock orientation. c These blending factors are based upon velocity jumps in order c to sense the direction of any compression, even if that c compression is not momentarily associated with a pressure c jump (as in the case of a shock reflecting from a rigid wall). c c c Now blend the 2. pre- and 2. post-shock values. c c do 3501 i = -mbdy+1,n+mbdy ppre(i) = factrx(i) * p0x(i) + factry(i) * p0y(i) + 1 factrz(i) * p0z(i) ppst(i) = factrx(i) * p1x(i) + factry(i) * p1y(i) + 1 factrz(i) * p1z(i) uxpre(i) = factrx(i) * ux0x(i) + factry(i) * ux0y(i) + 1 factrz(i) * ux0z(i) uxpst(i) = factrx(i) * ux1x(i) + factry(i) * ux1y(i) + 1 factrz(i) * ux1z(i) uypre(i) = factrx(i) * uy0x(i) + factry(i) * uy0y(i) + 1 factrz(i) * uy0z(i) 3501 continue do 3502 i = -mbdy+1,n+mbdy uypst(i) = factrx(i) * uy1x(i) + factry(i) * uy1y(i) + 1 factrz(i) * uy1z(i) uzpre(i) = factrx(i) * uz0x(i) + factry(i) * uz0y(i) + 1 factrz(i) * uz0z(i) uzpst(i) = factrx(i) * uz1x(i) + factry(i) * uz1y(i) + 1 factrz(i) * uz1z(i) rhopre(i) = factrx(i) * rho0x(i) + factry(i) * rho0y(i) + 1 factrz(i) * rho0z(i) rhopst(i) = factrx(i) * rho1x(i) + factry(i) * rho1y(i) + 1 factrz(i) * rho1z(i) 3502 continue c c do 3600 i = -mbdy+1,n+mbdy c c We may now proceed roughly as we do in 1. dimension. c First we turn off dffuse(i) outside shocks. c THE COMPUTATION OF divu BELOW ASSUMES A UNIFORM GRID. c pjump = ppst(i) - ppre(i) pjumpf = pjump / ppre(i) flag = shkjmp - pjumpf divu = ux(i+1) - ux(i-1) + uyzt(i) - uyzb(i) + uzzf(i) - uzzn(i) c c ***************************************************************** c OPPS. I THINK THAT DEFINING "NEAR" TO BE LOWER VALUES OF Z c ***************************************************************** c PRODUCES A LEFT-HANDED COORDINATE SYSTEM. c ***************************************************************** c shockd = max (divu, flag) c if ( shockd .gt. 0 ) shockd = 0. if ( shockd .lt. 0 ) shockd = 1. c dffuse(i) = shockd * dffuse(i) difuse(i) = dffuse(i) c c Now we control the time step for the special case of strong c shocks. Our formulae assume a single, straight shock front. c uxjump = uxpst(i) - uxpre(i) uyjump = uypst(i) - uypre(i) uzjump = uzpst(i) - uzpre(i) ujmpsq = uxjump * uxjump + uyjump * uyjump + uzjump * uzjump deenom = 1. / sqrt (ujmpsq + smlusq) w = pjump * deenom unitx = uxjump * deenom unity = uyjump * deenom unitz = uzjump * deenom c c unitx, unity, and unitz are the components of a unit vector point- c ing in a direction perpendicular to the shock front. In fact, c this unit vector points in the direction of shock propagation c relative to the fluid. c Note that the Lagrangian wave speed w above is always positive. c By constructing the velocities u0prp and u1prp using unitx, c unity, and unitz below, we transform all cases into coordinates c for which positive velocities are in the direction of shock motion c relative to the fluid. c u0prp = unitx * uxpre(i) + unity * uypre(i) + unitz * uzpre(i) c c not used?? This broke vectorization, removed as dead code c u1prp = unitx * uxpst(i) + unity * uypst(i) + unitz * uzpst(i) vpre = 1. / rhopre(i) weul = w * vpre c cepre = sqrt (gamma * ppre(i) * vpre) cepst = sqrt (gamma * ppst(i) / rhopst(i)) c cemax = max (cepre, cepst) weul = min (weul, cemax) wavspd = weul + u0prp c c WE WILL ASSUME THAT DY = DX. c Actually, I am not sure that any change need be made in the case c where dx and dy are not equal. c thyng = shockd * dtbydx(i) courwv = thyng * abs (wavspd) c c For the general equation of state code, cezbb and ceztt could c be made available to this routine as input arguments. c courno(i) = max (courno(i), courwv) 3600 continue c csuperscalar code c c c return end c c---------------------------------------------------------------------72 c subroutine dintrf (al, ar, da, a6, dal, dar, n, nbdy, noff) c c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** ***** c ***** dddd iii n n ttttt rrrr fffff ***** c ***** d d i nn n t r r f ***** c ***** d d i n n n t rrrr fff ***** c ***** d d i n n n t r r f ***** c ***** d d i n nn t r r f ***** c ***** dddd iii n n t r r f ***** c ***** ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c======================================================================= c dimension a6(-nbdy+1:n+nbdy+1), da(-nbdy+1:n+nbdy+1), 1 al(-nbdy+1:n+nbdy+1), ar(-nbdy+1:n+nbdy+1), 2 dal(-nbdy+1:n+nbdy+1), dar(-nbdy+1:n+nbdy+1) 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 March, 1993 ***** c ***** REVISED August, 1995 ***** c ***** Revised September, 1995 ***** c ***** ***** c ***** ***** c ***** COPYRIGHT Paul R. Woodward ***** c ***** March, 1994 ***** c ***** September, 1995 ***** c ***** ALL RIGHTS RESERVED ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c c c c This routine constructs the simplest possible monotone parabolae c within zones of the grid, using only the differences, dal, of c zone-averaged values of the variable. The array dar is used c for scratch storage space, and its values are not expected to c be of use to the calling program. dal and dar are, respectively, c differences of the zone averages centered at the left- and right- c hand interfaces of the zone. c When dal and dar have opposite signs, we are at an extremum, c and the interpolation parabola is set to a constant value a(i). c When both dal and dar are positive, the smaller of them is c always left unchanged, but the larger 1. may be reset to comply c with the monotonicity constraint that it be no larger than 4 times c the value of the smaller difference. c When both differences are negative, a similar constraint is c applied. c The interpolation parabola for each zone is constructed from these c constrained differences precisely as it would have been made from c dal and dar in the absence of any constraint. c c n is the number of real, physical zones in the grid. c nbdy is the number of fake zones at each end of the physical grid. c noff is the offset from each end of the full (real + fake) grid c which defines the strip in which meaningful values for c the ouput vectors al, ar, da, and a6 are to be c provided by this routine. c It is assumed that boundary conditions have been applied in the c calling program, so that all n+2*(nbdy-noff+1) zones of average c values are set properly. c c OUTPUT ARRAYS: c al = the interpolated value at the left-hand interface of the c zone, assuming a vanishing average value in the zone. c ar = the interpolated value at the right-hand interface of the c zone, assuming a vanishing average value in the zone. c da = ar - al c a6 = - 3 * (al + ar) c c Definitions of the above terms: c Let xtilde be a variable which assumes the values -1/2 and c 1/2 at the left- and right-hand interfaces of the zone of c interest. Also, we assume a uniform grid, so that the left-hand c interface of the zone to the left is at xtilde = -3/2, and the c right-hand interface of the sone on the right is at xtilde = 3/2. c Now, our interpolation parabola, which has average values in these c 3. zones equal to azl, a, and azr, the prescribed zone c averages, is given by the following formulae: c c a(xtilde) = am + da * xtilde - a6 * xtilde**2 c c al = a(-1/2) c ar = a(1/2) c c da = ar - al c abar = (al + ar) / 2 c a6 = 6 * (a - abar) c am = a + a6 / 12 c c If dal = a - azl c and dar = azr - a, c then: c dal = da + a6 c dar = da - a6 c Hence: c da = (dal + dar) / 2 c a6 = (dal - dar) / 2 c abar = a - (dal - dar) / 12 c al = a - (2*dal + dar) / 6 c ar = a + (2*dar + dal) / 6 c c Monotonicity is obtained by constraining the values of c dal and dar, and then using the formulae above to c obtain da, a6, abar, al, and ar. c If dal and dar have opposite signs, both are set to c 0.. Then, if both are positive, the larger of them is c constrained to be less than or equal to 4 times the smaller. c In the critical case, if dal were the smaller, then al c would just equal azl and (da/dxtilde) at the left-hand c interface would equal dal and the parabola would attain c its minimum value at xtilde = -5/6. c c mbdy = nbdy - noff c do 1000 i = -mbdy+1,n+mbdy+1 dar(i-1) = dal(i) 1000 continue c twelth = .25 / 3. c do 2000 i = -mbdy+1,n+mbdy c daal = dal(i) daar = dar(i) if ((daal * daar) .lt. 0.) then daal = 0 daar = 0 endif c if ((dal(i) * dar(i)) .lt. 0.) then c da(i) = 0. c a6(i) = 0. c al(i) = 0. c ar(i) = 0. c go to 2000 c else c daal = dal(i) c daar = dar(i) c endif c twodda = daal + daar c s = sign (1., twodda) s = 1. if ( twodda .lt. 0 ) s= -s daal = s * daal daar = s * daar daaal = s * min (daal, 4. * daar) daaar = s * min (daar, 4. * daal) da(i) = .5 * (daaal + daaar) abar = twelth * (daaar - daaal) a6(i) = - 6. * abar al(i) = abar - .5 * da(i) ar(i) = abar + .5 * da(i) 2000 continue c return end subroutine interf (a, al, ar, da, a6, dal, dar, n, nbdy, noff) c c c ************************************************** c ************************************************** c ************************************************** c ************************************************** c ***** ***** c ***** ***** c ***** iii n n ttttt eeee rrr fffff ***** c ***** i nn n t e r r f ***** c ***** i n n n t ee rrr fff ***** c ***** i n n n t e rr f ***** c ***** i n nn t e r r f ***** c ***** iii n n t eeee r r f ***** c ***** ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c======================================================================= c dimension a(-nbdy+1:n+nbdy+1), da(-nbdy+1:n+nbdy+1), 1 al(-nbdy+1:n+nbdy+1), ar(-nbdy+1:n+nbdy+1), 2 dal(-nbdy+1:n+nbdy+1), dar(-nbdy+1:n+nbdy+1), 3 a6(-nbdy+1:n+nbdy+1) 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 March, 1993 ***** c ***** REVISED August, 1995 ***** c ***** Revised September, 1995 ***** c ***** ***** c ***** ***** c ***** COPYRIGHT Paul R. Woodward ***** c ***** March, 1994 ***** c ***** September, 1995 ***** c ***** ALL RIGHTS RESERVED ***** c ***** ***** c ************************************************** c ************************************************** c ************************************************** c ************************************************** c c c c c This routine constructs the simplest possible monotone parabolae c within zones of the grid, using only the zone-averaged values, c a(i), of the variable. The 2. arrays dal and dar are used c for scratch storage space, and their values are not expected to c be of use to the calling program. They are, respectively, c differences of the zone averages centered at the left- and right- c hand interfaces of the zone. c When dal and dar have opposite signs, we are at an extremum, c and the interpolation parabola is set to a constant value a(i). c When both dal and dar are positive, the smaller of them is c always left unchanged, but the larger 1. may be reset to comply c with the monotonicity constraint that it be no larger than 4 times c the value of the smaller difference. c When both differences are negative, a similar constraint is c applied. c The interpolation parabola for each zone is constructed from these c constrained differences precisely as it would have been made from c dal and dar in the absence of any constraint. c c n is the number of real, physical zones in the grid. c nbdy is the number of fake zones at each end of the physical grid. c noff is the offset from each end of the full (real + fake) grid c which defines the strip in which meaningful values for c the ouput vectors al, ar, da, and a6 are to be c provided by this routine. c It is assumed that boundary conditions have been applied in the c calling program, so that all n+2*(nbdy-noff+1) zones of average c values are set properly. c c OUTPUT ARRAYS: c al = the interpolated value at the left-hand interface of the c zone, assuming a vanishing average value in the zone. c ar = the interpolated value at the right-hand interface of the c zone, assuming a vanishing average value in the zone. c da = ar - al c a6 = - 3 * (al + ar) c c Definitions of the above terms: c Let xtilde be a variable which assumes the values -1/2 and c 1/2 at the left- and right-hand interfaces of the zone of c interest. Also, we assume a uniform grid, so that the left-hand c interface of the zone to the left is at xtilde = -3/2, and the c right-hand interface of the sone on the right is at xtilde = 3/2. c Now, our interpolation parabola, which has average values in these c 3. zones equal to azl, a, and azr, the prescribed zone c averages, is given by the following formulae: c c a(xtilde) = am + da * xtilde - a6 * xtilde**2 c c al = a(-1/2) c ar = a(1/2) c c da = ar - al c abar = (al + ar) / 2 c a6 = 6 * (a - abar) c am = a + a6 / 12 c c If dal = a - azl c and dar = azr - a, c then: c dal = da + a6 c dar = da - a6 c Hence: c da = (dal + dar) / 2 c a6 = (dal - dar) / 2 c abar = a - (dal - dar) / 12 c al = a - (2*dal + dar) / 6 c ar = a + (2*dar + dal) / 6 c c Monotonicity is obtained by constraining the values of c dal and dar, and then using the formulae above to c obtain da, a6, abar, al, and ar. c If dal and dar have opposite signs, both are set to c 0.. Then, if both are positive, the larger of them is c constrained to be less than or equal to 4 times the smaller. c In the critical case, if dal were the smaller, then al c would just equal azl and (da/dxtilde) at the left-hand c interface would equal dal and the parabola would attain c its minimum value at xtilde = -5/6. c c mbdy = nbdy - noff c twelth = .25 / 3. c i = -mbdy+1 dal(i) = (a(i) - a(i-1)) dar(i-1) = dal(i) do 1000 i = -mbdy+1,n+mbdy dal(i+1) = (a(i+1) - a(i)) dar(i) = dal(i+1) 1000 continue c do 2000 i = -mbdy+1,n+mbdy c daal = dal(i) daar = dar(i) if ( daal * daar .lt. 0 ) then daal = 0 daar = 0 endif c if ((dal(i) * dar(i)) .lt. 0.) then c da(i) = 0. c a6(i) = 0. c al(i) = a(i) c ar(i) = a(i) c go to 2000 c else c daal = dal(i) c daar = dar(i) c endif c twodda = daal + daar c s = sign (1., twodda) s = 1. if ( twodda .lt. 0 ) s= -s daal = s * daal daar = s * daar daaal = s * min (daal, 4. * daar) daaar = s * min (daar, 4. * daal) da(i) = .5 * (daaal + daaar) abar = a(i) + twelth * (daaar - daaal) a6(i) = 6. * (a(i) - abar) al(i) = abar - .5 * da(i) ar(i) = abar + .5 * da(i) 2000 continue c return end c c csuncode( c function cvmgm (a, b, f) cvmgm = b if (f .lt. 0.) cvmgm = a return end function cvmgz (a, b, f) cvmgz = b if (f .eq. 0.) cvmgz = a return end c csuncode) c