PROGRAM MAP_2 c c Givn two uniform grids defined by the locations of the lefthand c zone interfaces XL and XNUL. The grids are offset by an amount c DX, that is XNUL(I) = XL(I) + DX. c This program creates piecewise parabolic interpolations c from zone averaged data, defined on the XL grid. c Monotonicity constraints are applied. For those zones c considered to be in a discontinuity, an interpolation more c appropriate for discontinuous functions is used, producing c a better representation of unresolved shocks. c Regions of over lap between the two grid are determined and c advective fluxes found by integrating the PPM interpolations c over these regions. The zone averages on the grid XNUL are found c by differeing these fluxes. c c Note, the interpolations require 3 boundary zones at each end, and c the integrations will require an additional zone. c PARAMETER ( NBDY = 4 ) ALLOCATABLE :: XL(:), A_AVG(:), AL(:), AR(:), DA(:), A6(:), ^ SMALLDA(:), UNSMOOTH(:), DISCONT(:), DX(:), DA_AVL(:), ^ DVOLL(:), XF(:), XF1(:), XLNU(:), A_AVGNU(:) c OPEN (unit = 11, file ='INPUT.AVG') c The data consists of N+2 lines. The first line c contains the integer N, the number to zone averages in the data set. c This is followed by N lines, each line having two numbers, a value for c the location of the left edge of the zone XL, and the zone average for c that zone A_AVG. The lines are assumed to be in sequential order according c to XL. An additional line has only one value, the final left zone edge. A c grid of N zones requires N+1 left zone edge locations. An appropriate data c set can be created from the program CREATE_ZONE_AVERAGES. The XL values c are used for plotting, not in the interpolation process. c Read in the number of zones in the file READ (11,*) N c c Next, allocate the appropriate space ALLOCATE (XL(1-NBDY:N+1+NBDY)) ALLOCATE (AL(1-NBDY:N+1+NBDY)) ALLOCATE (AR(1-NBDY:N+1+NBDY)) ALLOCATE (DA(1-NBDY:N+NBDY)) ALLOCATE (A6(1-NBDY:N+NBDY)) ALLOCATE (A_AVG(1-NBDY:N+NBDY)) ALLOCATE (DA_AVL(1-NBDY:N+NBDY)) ALLOCATE (XLNU(1-NBDY:N+1+NBDY)) ALLOCATE (A_AVGNU(1-NBDY:N+NBDY)) ALLOCATE (SMALLDA(1-NBDY:N+NBDY)) ALLOCATE (UNSMOOTH(1-NBDY:N+NBDY)) ALLOCATE (DISCONT(1-NBDY:N+NBDY)) ALLOCATE (DVOLL(1-NBDY:N+NBDY)) ALLOCATE (XF(1-NBDY:N+NBDY)) ALLOCATE (XF1(1-NBDY:N+NBDY)) c The following arrays will hold various geometrical factors c same grid: ALLOCATE (DX(1-NBDY:N+NBDY)) DO I=1,N READ (11,*) XL(I),A_AVG(I) ENDDO READ (11,*) XL(N+1) c c c Query the user for desired boundary conditions c PRINT *,' Please enter 0 for reflecting boundary conditions' PRINT *,' or 1 for periodic boundary conditions' READ *, IBDY_TYPE IF ( IBDY_TYPE .EQ. 0 ) THEN PRINT *,' Selecting reflecting boundary conditions' ELSE IF ( IBDY_TYPE .EQ. 1 ) THEN PRINT *,' Selecting periodic boundary conditions' ELSE PRINT *,' Unknown boundary type:',IBDY_TYPE call abort() ENDIF ENDIF c c Initialize the fake zones: CALL BOUNDARY( XL(1-NBDY), A_AVG(1-NBDY), IBDY_TYPE, N,NBDY ) c c c While a uniform grid is assumed, it is convient to calculate c the individual zone widths. DO I=1-NBDY,N+NBDY DX(I) = XL(I+1) - XL(I) ENDDO c c Now determine the offset between the two grids. For stability, we require c that the offset be no greater than one zone. That is c XL(I-1) <= XLNU(I) <= XL(I+1), when express as a fraction COURNO c c COURNO = ABS(XL(I) - XLNU(I))/DX c this requirement becores c c 0 <= COURNO(I) <= 1.0, c c A familiar Courant condition or stability condition. c c 1 PRINT *,' Enter the fractional offset (between -1.0 and 1.0)' PRINT *,' between the two grids: DXF. ' READ *,DXF IF( DXF .LT. -1.0 .OR. DXF .GT. 1.0 ) THEN PRINT *,' You entered',DXF GOTO 1 ENDIF c c Determine the locations of the lefthand zone interface of the new grid c (The first and last interface locations will not be needed) c DO I=2-NBDY,N+NBDY IF (DXF .LE. 0.0) THEN XLNU(I) = XL(I) + DXF * DX(I-1) ELSE XLNU(I) = XL(I) + DXF * DX(I) ENDIF ENDDO c Fill an array, SMALLDA, which contains what is to be considered a trivial difference of A_AVG. c In this Program, A_AVG is not positive definite, so a small value is appropriate. DO I=1-NBDY,N+NBDY SMALLDA(I) = 0.0005 ENDDO c Interpolate zone interface values and determine ppm parabola coefficients c CALL PPM98_INTERP0_DISCONTINUITY ^ ( SMALLDA(1-NBDY), A_AVG(1-NBDY), AL(1-NBDY), ^ AR(1-NBDY), DA(1-NBDY), A6(1-NBDY), ^ UNSMOOTH(1-NBDY), DISCONT(1-NBDY), N, NBDY ) c AL(I), AR(I), DA(I), A6(I), are defined for 4-NBDY<= I <= N+NBDY-3 FORTHD = 4.0 / 3.0 nnbdy = NBDY-3 DO I=1-nnbdy,N+nnbdy DVOLL(i) = XL(I) - XLNU(I) IF (DVOLL(I) .GE. 0.0) THEN c The data will come from the c zone (i-1) at the left, c XF(I) = DVOLL(I)/DX(i-1) c ELSE c c The data will come from zone (i) XF(I) = -DVOLL(I)/DX(I) c ENDIF XF1(I) = 1.0 - XF(I) ENDDO mbdy = 1 CALL PPM98_FLUX(AL(1-NNBDY), A_AVG(1-NNBDY), AR(1-NNBDY), 6 DA(1-NNBDY), DVOLL(1-NNBDY), XF(1-NNBDY), ^ XF1(1-NNBDY), DA_AVL(1-NNBDY), N, MBDY) c c Now use the fluxes DA_AVL to find the new zone averages c DO I=5-NBDY,N+NBDY-4 DXNU = XLNU(I+1)-XLNU(I) A_AVGNU(I) = ( A_AVG(I)*DX(I) + DA_AVL(I)-DA_AVL(i+1))/DXNU ENDDO c c Display old and new zone averages PRINT *,' I XMID A_AVG XMIDNU A_AVGNU' DO I=5-NBDY,N+NBDY-4 XMID = 0.5*(XL(I+1)+XL(I)) XMIDNU = 0.5*(XLNU(I+1)+XLNU(I)) WRITE(6,100) I, XMID, A_AVG(I), XMIDNU, A_AVGNU(I) 100 FORMAT(3x,i4,1p,4(1x,e12.5)) ENDDO DEALLOCATE (XL, A_AVG, AL, AR, DA, A6, SMALLDA, UNSMOOTH, DISCONT) DEALLOCATE (DX, A_AVGNU, XLNU, DA_AVL, DVOLL, XF, XF1) END ! PROGRAM MAP_2