PROGRAM ADVECT_1 c c This program developed from the INTERPOLATE_4 example. c c This program creates piecewise parabolic interpolations c from zone averaged data, without assuming a uniform 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. A second c non-uniform grid (grid B) offset (by less than one zone) from c the initial grid (grid A) is created. Appropriate fluxes c are computed which will allow the determination of the zone c averages on grid B using the zone averages on the grid A. c These fluxes are found by integrating the piecewise c parabolic interpolations (determined on grid A) over the c regions of overlap between the two grids. c c The number of boundary zones is: 3 for interpolation and 1 for c integration = 4. To determine the lefthand zone interface flux c for zone I, the interpolation in zone I-1 may be required, c depending upon zone motion. Likewise, to determine the c righthand zone interface flux for zone I, the interpolation in c zone I+1 may be required. c PARAMETER ( NBDY = 4 ) ALLOCATABLE :: XL(:), A_AVG(:), AL(:), AR(:), DA(:), A6(:), ^ SMALLDA(:), UNSMOOTH(:), DISCONT(:), DX(:), DXINV(:), ^ XLNU(:), A_AVGNU(:), DVOLL(:), DVOLFL(:), DVOLFL1(:), ^ DAL(:) ALLOCATABLE :: DALFAC(:), DARFAC(:), FACDA(:), FCDAZL(:), ^ FACDAL(:), FA6DAL(:), FA6DAR(:), FRRDAL(:), FRRDAR(:), ^ FLLDAL(:), FLLDAR(:) 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 (SMALLDA(1-NBDY:N+NBDY)) ALLOCATE (UNSMOOTH(1-NBDY:N+NBDY)) ALLOCATE (DISCONT(1-NBDY:N+NBDY)) ALLOCATE (DX(1-NBDY:N+NBDY)) ALLOCATE (DXINV(1-NBDY:N+NBDY)) ALLOCATE (XLNU(1-NBDY:N+NBDY)) ALLOCATE (A_AVGNU(1-NBDY:N+NBDY)) ALLOCATE (DVOLL(1-NBDY:N+NBDY)) ALLOCATE (DVOLFL(1-NBDY:N+NBDY)) ALLOCATE (DVOLFL1(1-NBDY:N+NBDY)) ALLOCATE (DAL(1-NBDY:N+NBDY)) c The following arrays will hold various geometrical factors which can be c precomputed for efficiency when several quantities are interpolated on the c same grid: ALLOCATE (DALFAC(1-NBDY:N+NBDY)) ALLOCATE (DARFAC(1-NBDY:N+NBDY)) ALLOCATE (FACDA(1-NBDY:N+NBDY)) ALLOCATE (FCDAZL(1-NBDY:N+NBDY)) ALLOCATE (FACDAL(1-NBDY:N+NBDY)) ALLOCATE (FA6DAL(1-NBDY:N+NBDY)) ALLOCATE (FA6DAR(1-NBDY:N+NBDY)) ALLOCATE (FRRDAL(1-NBDY:N+NBDY)) ALLOCATE (FRRDAR(1-NBDY:N+NBDY)) ALLOCATE (FLLDAL(1-NBDY:N+NBDY)) ALLOCATE (FLLDAR(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 A uniform grid is NOT assumed. Therefore, various c geometrical factors must be precomputed by a call to c _INTERFACE_FACTORS, input to which consists of DX(I), c the individual zone widths. DO I=1-NBDY,N+NBDY DX(I) = XL(I+1) - XL(I) DXINV(I) = 1.0 / DX(I) ENDDO CALL PPM98_INTERPOLATION_FACTORS ( DX(1-NBDY), DALFAC(1-NBDY), ^ DARFAC(1-NBDY), FACDA(1-NBDY), FCDAZL(1-NBDY), ^ FACDAL(1-NBDY), FA6DAL(1-NBDY), FA6DAR(1-NBDY), ^ FRRDAL(1-NBDY), FRRDAR(1-NBDY), FLLDAL(1-NBDY), ^ FLLDAR(1-NBDY), N, NBDY) c 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 c Interpolate zone interface values and determine ppm parabola coefficients c CALL PPM98_INTERP_DISCONTINUITY ( DALFAC(1-NBDY), DARFAC(1-NBDY), ^ FACDA(1-NBDY), FCDAZL(1-NBDY), FACDAL(1-NBDY), ^ FA6DAL(1-NBDY), FA6DAR(1-NBDY), FLLDAL(1-NBDY), ^ FLLDAR(1-NBDY), FRRDAL(1-NBDY), FRRDAR(1-NBDY), ^ 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 c 1 PRINT *,' Graciously enter the fractional offset (between -1.0 and 1.0)' PRINT *,' between the old and new grids: DXF. ' READ *,DXF IF( DXF .LT. -1.0 .OR. DXF .GT. 1.0 ) THEN PRINT *,' You entered',DXF, ' Do please try again.' GOTO 1 ENDIF PRINT *,' The grids will be separated by ',DXF,' times ' PRINT *,' the smallest zone width. Thank you.' c c Find the minimum zone width: c DXMIN = DX(1-NBDY) DO I=2-NBDY,N+NBDY DXMIN = AMIN1( DXMIN, DX(I)) ENDDO c c Determine the new grid c DO i=2-NBDY,N+NBDY XLNU(I) = XL(I) + DXF * DXMIN ENDDO c c While the following operations are trivial given the constant c offset, they are explicitly computed here to show how the c general case would proceed c c Find advected volumes at left hand zone interfaces, c also find advected volume fractions at left hand zone interfaces. DO i=2-NBDY,N+NBDY c c Advected volume: DVOLL(I) = XLNU(I) - XL(I) c Advected volume fraction: IF ( DVOLL(I) .LE. 0.0 ) THEN DVOLFL(I) = - DVOLL(I) * DXINV(I-1) ELSE DVOLFL(I) = DVOLL(I) * DXINV(I) ENDIF DVOLFL1(I) = 1.0 - DVOLFL(I) ENDDO c c Find the fluxes of lefthand zone interfaces fluxes of A (DAL) which will c map zone averages of A from the grid XL to the grid XLNU: c c Modify offsets and boundary zone count because 3 fake zones have c been lost to the interpolation. c NNBDY = NBDY-3 MBDY = NBDY-3 c CALL PPM98_FLUX( AL(1-NNBDY), A_AVG(1-NNBDY), AR(1-NNBDY), DA(1-NNBDY), ^ DVOLL(1-NNBDY), DVOLFL(1-NNBDY), DVOLFL1(1-NNBDY), DAL(1-NNBDY), ^ N, MBDY ) c c Now determine the new zone averages c DO I=1,N DXNEW = XLNU(I+1)-XLNU(I) print *,' dx,dxnu=',dx(i),dxnew,dxnew*dxinv(i) A_AVGNU(I) = ( A_AVG(I) * DX(I) + DAL(I) - DAL(I+1) ) / DXNEW ENDDO c c Display the old and new zone averages PRINT *,' I XMID A_AVG XMIDNU A_AVGNU' DO I=1,N XMID = 0.5*(XL(I)+XL(I+1)) XMIDNU = 0.5*(XLNU(I)+XLNU(I+1)) WRITE(6,100) I, XMID, A_AVG(I), XMIDNU, A_AVGNU(I) ENDDO 100 FORMAT( I4,1P, 4(1X,e12.5)) DEALLOCATE ( XL, A_AVG, AL, AR, DA, A6, ^ SMALLDA, UNSMOOTH, DISCONT, DX, DXINV, ^ XLNU, A_AVGNU, DVOLL, DVOLFL, DVOLFL1, ^ DAL ) END ! PROGRAM ADVECT_1