Page 1 (printed 3/30/99) DDDDOOOO____PPPPPPPPMMMMDDDDEEEE0000____1111DDDDCCCC____GGGGAAAAMMMMMMMMAAAA((((JJJJaaaannnn 1111999999999999))))UUUUNNNNIIIIXXXX SSSSyyyysssstttteeeemmmm VVVVDDDDOOOO____PPPPPPPPMMMMDDDDEEEE0000____1111DDDDCCCC____GGGGAAAAMMMMMMMMAAAA((((JJJJaaaannnn 1111999999999999)))) NNNNAAAAMMMMEEEE ddddoooo____ppppppppmmmmddddeeee0000____1111ddddcccc____ggggaaaammmmmmmmaaaa - Given 1d strips of zone averages of density, pressure, and velocity defined on a uniform grid, new zone averages of these quantities are returned by advancing the Euler equations forward one timestep. FFFFOOOORRRRTTTTRRRRAAAANNNN SSSSPPPPEEEECCCCIIIIFFFFIIIICCCCAAAATTTTIIIIOOOONNNN ssssuuuubbbbrrrroooouuuuttttiiiinnnneeee ddddoooo____ppppppppmmmmddddeeee0000____1111ddddcccc____ggggaaaammmmmmmmaaaa(((( XXXXLLLL,,,, RRRRHHHHOOOO,,,, PPPP,,,, UUUU,,,, RRRRHHHHOOOONNNNUUUU,,,, PPPPNNNNUUUU,,,, UUUUNNNNUUUU,,,, EEEEOOOOSSSSGGGGAAAAMMMM,,,, DDDDTTTT,,,, SSSSMMMMLLLLRRRRHHHHOOOO,,,, SSSSMMMMAAAALLLLLLLLPPPP,,,, SSSSMMMMAAAALLLLLLLLUUUU,,,, SSSSMMMMAAAALLLLLLLLEEEE,,,, CCCCOOOOUUUURRRRMMMMXXXX,,,, NNNNZZZZOOOONNNNEEEESSSS,,,, NNNNBBBBDDDDYYYY )))) rrrreeeeaaaallll****4444 XXXXLLLL((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY++++1111)))) rrrreeeeaaaallll****4444 RRRRHHHHOOOO((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY)))) rrrreeeeaaaallll****4444 PPPP((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY)))) rrrreeeeaaaallll****4444 UUUU((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY)))) rrrreeeeaaaallll****4444 RRRRHHHHOOOONNNNUUUU((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY)))) rrrreeeeaaaallll****4444 PPPPNNNNUUUU((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY)))) rrrreeeeaaaallll****4444 UUUUNNNNUUUU((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY)))) rrrreeeeaaaallll****4444 EEEEOOOOSSSSGGGGAAAAMMMM,,,, DDDDTTTT,,,, SSSSMMMMLLLLRRRRHHHHOOOO,,,, SSSSMMMMAAAALLLLLLLLPPPP.... SSSSMMMMAAAALLLLLLLLUUUU,,,, SSSSMMMMAAAALLLLLLLLEEEE,,,, rrrreeeeaaaallll****4444 CCCCOOOOUUUURRRRMMMMXXXX iiiinnnntttteeeeggggeeeerrrr****4444 NNNNZZZZOOOONNNNEEEESSSS,,,, NNNNBBBBDDDDYYYY ssssuuuubbbbrrrroooouuuuttttiiiinnnneeee dddd____ddddoooo____ppppppppmmmmddddeeee0000____1111ddddcccc____ggggaaaammmmmmmmaaaa(((( XXXXLLLL,,,, RRRRHHHHOOOO,,,, PPPP,,,, UUUU,,,, RRRRHHHHOOOONNNNUUUU,,,, PPPPNNNNUUUU,,,, UUUUNNNNUUUU,,,, EEEEOOOOSSSSGGGGAAAAMMMM,,,, DDDDTTTT,,,, SSSSMMMMLLLLRRRRHHHHOOOO,,,, SSSSMMMMAAAALLLLLLLLPPPP,,,, SSSSMMMMAAAALLLLLLLLUUUU,,,, SSSSMMMMAAAALLLLLLLLEEEE,,,, CCCCOOOOUUUURRRRMMMMXXXX,,,, NNNNZZZZOOOONNNNEEEESSSS,,,, NNNNBBBBDDDDYYYY )))) rrrreeeeaaaallll****8888 XXXXLLLL((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY++++1111)))) rrrreeeeaaaallll****8888 RRRRHHHHOOOO((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY)))) rrrreeeeaaaallll****8888 PPPP((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY)))) rrrreeeeaaaallll****8888 UUUU((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY)))) rrrreeeeaaaallll****8888 RRRRHHHHOOOONNNNUUUU((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY)))) rrrreeeeaaaallll****8888 PPPPNNNNUUUU((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY)))) rrrreeeeaaaallll****8888 UUUUNNNNUUUU((((1111----NNNNBBBBDDDDYYYY::::NNNNZZZZOOOONNNNEEEESSSS++++NNNNBBBBDDDDYYYY)))) rrrreeeeaaaallll****8888 EEEEOOOOSSSSGGGGAAAAMMMM,,,, DDDDTTTT,,,, SSSSMMMMLLLLRRRRHHHHOOOO,,,, SSSSMMMMAAAALLLLLLLLPPPP.... SSSSMMMMAAAALLLLLLLLUUUU,,,, SSSSMMMMAAAALLLLLLLLEEEE,,,, rrrreeeeaaaallll****8888 CCCCOOOOUUUURRRRMMMMXXXX iiiinnnntttteeeeggggeeeerrrr****4444 NNNNZZZZOOOONNNNEEEESSSS,,,, NNNNBBBBDDDDYYYY These two routines are functionally the same. They differ only in the type declarations of their (real) arguments (that is, single or double precision.) AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS XXXXLLLL [input] array of the locations of the lefthand zone Page 1 (printed 3/30/99) DDDDOOOO____PPPPPPPPMMMMDDDDEEEE0000____1111DDDDCCCC____GGGGAAAAMMMMMMMMAAAA((((JJJJaaaannnn 1111999999999999))))UUUUNNNNIIIIXXXX SSSSyyyysssstttteeeemmmm VVVVDDDDOOOO____PPPPPPPPMMMMDDDDEEEE0000____1111DDDDCCCC____GGGGAAAAMMMMMMMMAAAA((((JJJJaaaannnn 1111999999999999)))) interfaces. This defines the grid. Values on the range [1-NBDY, NZONES+NBDY+1] are required RRRRHHHHOOOO [input] array of zone averages of density. Values on the range [1-NBDY, NZONES+NBDY] are required PPPP [input] array of zone averages of pressure. Values on the range [1-NBDY, NZONES+NBDY] are required UUUU [input] array of zone averages of velocity. Values on the range [1-NBDY, NZONES+NBDY] are required RRRRHHHHOOOONNNNUUUU [output] array of updated zone averages of density. Values on the range [1, NZONES] are computed PPPPNNNNUUUU [output] array of updated zone averages of pressure. Values on the range [1, NZONES] are computed UUUUNNNNUUUU [output] array of updated zone averages of velocity. Values on the range [1, NZONES] are computed EEEEOOOOSSSSGGGGAAAAMMMM [input] the ratio of specific heats. A gamma-law equation of state is used in which p(i)=(gamma- 1)*rho(i)*ei(i), with ei being the specific internal energy. DDDDTTTT [input] the timestep. SSSSMMMMLLLLRRRRHHHHOOOO [input] This indicates what should be considered a trivial values of density. It is used to reduce the effects of roundoff errors. Reasonable values are 1.0e-06 for single precision and 1.0e-08 for double precision. Page 2 (printed 3/30/99) DDDDOOOO____PPPPPPPPMMMMDDDDEEEE0000____1111DDDDCCCC____GGGGAAAAMMMMMMMMAAAA((((JJJJaaaannnn 1111999999999999))))UUUUNNNNIIIIXXXX SSSSyyyysssstttteeeemmmm VVVVDDDDOOOO____PPPPPPPPMMMMDDDDEEEE0000____1111DDDDCCCC____GGGGAAAAMMMMMMMMAAAA((((JJJJaaaannnn 1111999999999999)))) SSSSMMMMAAAALLLLLLLLPPPP [input] This indicates what should be considered a trivial values of pressure. It is used to reduce the effects of roundoff errors. Reasonable values are 1.0e-06 for single precision and 1.0e-08 for double precision. SSSSMMMMAAAALLLLLLLLUUUU [input] This indicates what should be considered a trivial values of velocity. It is used to reduce the effects of roundoff errors. Reasonable values are 1.0e-06 for single precision and 1.0e-08 for double precision. SSSSMMMMAAAALLLLLLLLEEEE [input] This indicates what should be considered a trivial values of energy. It is used to reduce the effects of roundoff errors. Reasonable values are 1.0e-06 for single precision and 1.0e-08 for double precision. CCCCOOOOUUUURRRRMMMMXXXX [input] The maximum Courant number encountered during this call. The original value of COURMX is lost during the computation. The Courant number is a stability condition and should be less than one. NNNNZZZZOOOONNNNEEEESSSS [input] the array dimension, or number of zones to process, not counting boundary or ``fake'' zones. NNNNBBBBDDDDYYYY [input] the number of boundary of ``fake'' zones which have been set at each end of the 1d arrays to implement the users choice of boundary conditions. This must be atleast 5. DDDDEEEESSSSCCCCRRRRIIIIPPPPTTTTIIIIOOOONNNN These routines are examples of one-dimensional direct Eulerian hydrodynamic subroutines using elements from the PPMLIB library. First, piecewise Parabolic interpolations are constructed for pressure and velocity. Where the input data distribution is not well represented by parabolae or poorly resolved, (that is, the data are ``unsmooth''), monotonicity constraints are applied. Piecewise Parabolic Page 3 (printed 3/30/99) DDDDOOOO____PPPPPPPPMMMMDDDDEEEE0000____1111DDDDCCCC____GGGGAAAAMMMMMMMMAAAA((((JJJJaaaannnn 1111999999999999))))UUUUNNNNIIIIXXXX SSSSyyyysssstttteeeemmmm VVVVDDDDOOOO____PPPPPPPPMMMMDDDDEEEE0000____1111DDDDCCCC____GGGGAAAAMMMMMMMMAAAA((((JJJJaaaannnn 1111999999999999)))) interpolations are also constructed for the density. Contact discontinuity detection is performed and where it is determined likely that a contact discontinuity exists, the density interpolations are modified to better represent discontinuities, rather than a smooth distribution. Monotonicity constrains are applied where the density data are deemed unsmooth. Second, appropriate left and right states are found for a Riemann problem at each zone interface. These state are calculated by integrating the above interpolations over the left and right domains of dependence corresponding to each lefthand zone interface. This results in a miniature Riemann problem, or discontinuity at each zone interface. Next, a Riemann solver is used to resolve these resultant discontinuities into time-averaged fluxes at the zone interfaces. These fluxes are differenced conservatively to provide new zone-averaged values of density, pressure (or energy), and velocity. Finally, a diffusion step is applied. Shock detection is performed, so that the additional diffusion is applied only in those zone which require it. This includes strong shocks and certain pathological case such as a shock moving slowly through the numerical grid. In addition to the updated zone averages of density, pressure, and velocity, an estimate of the Courant number is also returned to the calling program. In order for the method to be stable, the various Domains of Dependence for each zone interface may not exceed one zone in width. Nor may any wave travel more than one zone width during the timestep. To be stable, the Courant number should be less than 1.0. A gamma law equation of state is assumed: p = ( gamma -1 ) * rho * ei where p is the pressure, rho is the density, ei is the internal energy, and gamma is the ratio of specific heats. Gamma is assumed to be constant and uniform. HHHHIIIINNNNTTTTSSSS +o The required problem boundary conditions may be enforced by setting appropriate zone-averaged values of the fake zones. Page 4 (printed 3/30/99) DDDDOOOO____PPPPPPPPMMMMDDDDEEEE0000____1111DDDDCCCC____GGGGAAAAMMMMMMMMAAAA((((JJJJaaaannnn 1111999999999999))))UUUUNNNNIIIIXXXX SSSSyyyysssstttteeeemmmm VVVVDDDDOOOO____PPPPPPPPMMMMDDDDEEEE0000____1111DDDDCCCC____GGGGAAAAMMMMMMMMAAAA((((JJJJaaaannnn 1111999999999999)))) +o The timestep should be adjusted so that the Courant number (COURMX) remains less than 1.0. Because the actual Courant number encountered during a call is not known until after the call is completed, a common procedure is to adjust the new timestep so that the Courant number from the just completed call is <= 0.8. +o For an example of calls to the general Equation of state routines, see [[[[dddd____]]]]ddddoooo____ppppppppmmmmddddeeee0000____1111ddddcccc. +o For an example of an non uniform grid see [[[[dddd____]]]]ddddoooo____ppppppppmmmmddddeeee____1111ddddcccc____ggggaaaammmmmmmmaaaa. +o For examples of higher dimensions see [[[[dddd____]]]]ddddoooo____ppppppppmmmmddddeeee0000____2222ddddcccc____ggggaaaammmmmmmmaaaa or [[[[dddd____]]]]ddddoooo____ppppppppmmmmddddeeee0000____3333ddddcccc____ggggaaaammmmmmmmaaaa . BBBBUUUUGGGGSSSS +o SSSSEEEEEEEE AAAALLLLSSSSOOOO DO_PPMDE0_1DC, DO_PPMDE_1DC_GAMMA, DO_PPMDE0_2DC_GAMMA, DO_PPMDE0_3DC_GAMMA, FFFFIIIILLLLEEEESSSS ppm98_dfns.h ppm98_bdys.h BBBBUUUUGGGG RRRREEEEPPPPOOOORRRRTTTTSSSS TTTTOOOO ppmlib@sapphire.lcse.umn.edu AAAAUUUUTTTTHHHHOOOORRRRSSSS PPPPaaaauuuullll RRRR.... WWWWooooooooddddwwwwaaaarrrrdddd paul@lcse.umn.edu BBBB.... KKKKeeeevvvviiiinnnn EEEEddddggggaaaarrrr bke@lcse.umn.edu Department of Astronomy University of Minnesota Minneapolis, MN 55455 USA Page 5 (printed 3/30/99)