!=================================================================== ! multi1D.f90 - Solve the 1D Euler equations using Harten's ! multiresolution algorithm ! ! Ref.: A. Harten "Multiresolution Algorithms for the Numerical ! Solution of Hyperbolic Conservation Laws" ! Comm. on Pure and Applied Math. 48:1305-1342 ! (1995) ! ! Sorin Mitran Jan. 25, 2000 !=================================================================== PROGRAM Multi1D INTEGER p REAL, ALLOCATABLE, DIMENSION(:,:) :: vn,eps CALL CheckEnDeCode CALL ReadData(N,nbc,mcomp,L) ALLOCATE(vn(1-nbc:N+nbc,mcomp),eps(L,mcomp),STAT=ierror) IF (ierror .NE. 0) THEN PRINT 1001,ierror 1001 FORMAT('Error ',i6,' in allocating vn, eps arrays.') STOP END IF CALL Init(N,nbc,mcomp,vn,delt,delx,tfinal,CFL,nFlux,Kflux,nPlot,p,L,eps) CALL TimeStep(N,nbc,mcomp,vn,delt,delx,tfinal,CFL,nFlux,Kflux,nPlot,p,L,eps) CALL Fini DEALLOCATE(vn,eps) END PROGRAM Multi1D ! Check the encoding and decoding procedures with a simple sine function ! or a step function SUBROUTINE CheckEnDecode PARAMETER (L=4,N=2**L,sInterp=1) REAL v(0:N),vsave(0:N) INTEGER s COMMON /Interp/s delx=2.0/N x=-1.0 pi=4.0*atan(1.0) s=sInterp DO j=0,N v(j)=sin(pi*x) x=x+delx END DO ! Remove(Place) comments below for test of en/de-coding with a step(sine) function DO j=0,N/2+1 v(j)=1 x=x+delx END DO DO j=N/2+2,N v(j)=0 x=x+delx END DO vsave=v DO level=0,L v=vsave errmax=0.0 CALL PointEncode(N,v,level) CALL PointDecode(N,v,level) DO j=0,N err=ABS(v(j)-vsave(j)) IF (err > errmax) errmax=err END DO PRINT 1001,level,errmax 1001 FORMAT('Maximum error in point en/de-coding procedure to level ',i2, & ' is ',e10.3) END DO DO level=0,L v=vsave errmax=0.0 CALL CellEncode(N,v(1),level) CALL CellDecode(N,v(1),level) DO j=1,N err=ABS(v(j)-vsave(j)) IF (err > errmax) errmax=err END DO PRINT 1002,level,errmax 1002 FORMAT('Maximum error in cell en/de-coding procedure to level ',i2, & ' is ',e10.3) END DO END SUBROUTINE CheckEnDecode SUBROUTINE ReadData(N,nbc,mcomp,L) OPEN(UNIT=1,FILE='multi.dat',STATUS='OLD') READ(1,*)N READ(1,*)nbc READ(1,*)mcomp READ(1,*)L CLOSE(UNIT=1) END SUBROUTINE ReadData SUBROUTINE Init(N,nbc,mcomp,vn,delt,delx,tfinal,CFL,nFlux,Kflux,nPlot,p,L,eps) INTEGER p REAL vn(1-nbc:N+nbc,mcomp),eps(L,mcomp) REAL, ALLOCATABLE, DIMENSION(:) :: vleft,vright ! Read problem data OPEN(UNIT=1,FILE='multi.dat',STATUS='OLD') READ(1,*)N READ(1,*)nbc READ(1,*)mcomp READ(1,*)L READ(1,*)delt READ(1,*)delx READ(1,*)tfinal READ(1,*)CFL READ(1,*)nFlux DO i=1,3 READ(1,*) ! Skip comment lines for flux codes END DO READ(1,*)Kflux READ(1,*)nPlot READ(1,*)p ! Read left and right states of shock tube problem ALLOCATE(vleft(mcomp),vright(mcomp),STAT=ierror) IF (ierror .NE. 0) THEN PRINT 1001,ierror 1001 FORMAT('Error ',i6,' in allocating vn, eps arrays.') STOP END IF READ(1,*)jface READ(1,*)(vleft(m),m=1,mcomp) READ(1,*)(vright(m),m=1,mcomp) ! Truncation criteria on each level, for each component DO k=1,L READ(1,*)(eps(k,m),m=1,mcomp) END DO CLOSE(UNIT=1) ! Initial data DO m=1,mcomp vn(1-nbc:jface,m)=vleft(m) vn(jface+1:N+nbc,m)=vright(m) END DO OPEN(UNIT=2,FILE='multi.out',STATUS='REPLACE') ! Plot file WRITE(2,*)N,mcomp END SUBROUTINE Init ! The specific interpolation procedure used for function ! point values is encapsulated in this function REAL FUNCTION PointInterp(fk,j,N) REAL fk(0:N) INTEGER j,N REAL beta(2,2) INTEGER s ! The interpolation order COMMON /Interp/s DATA beta/ 0.50000000e0, 0.000000e0, & 0.56250000e0,-0.062500e0 / sum=0.0 DO l=1,s sum=sum+beta(l,s)*(fk(j+l-1)+fk(j-l)) END DO PointInterp=sum END FUNCTION PointInterp ! The specific interpolation procedure used for function ! primitive values is encapsulated in this function ! This should be consistent with the point interpolation procedure REAL FUNCTION PrimitiveInterp(fk,j,N) REAL fk(N) INTEGER j,N REAL gamma(2,2) INTEGER s ! The point interpolation order COMMON /Interp/s DATA gamma/ -0.12500000e0, 0.00000000e0, & -0.17187500e0, 0.23438000e0 / sum=0.0 DO l=1,s IF ((j+l) <= N) THEN fp=fk(j+l) ELSE fp=0.0 END IF IF ((j-l) >= 1) THEN fm=fk(j-l) ELSE fm=0.0 END IF sum=sum+gamma(l,s)*(fp-fm) END DO PrimitiveInterp=sum END FUNCTION PrimitiveInterp ! Function to encode function cell average values SUBROUTINE CellEncode(N,v,L) INTEGER N,L REAL v(N) ! Input: N - number of cells, must be a multiple of 2**L ! v - cell average values of function to be encoded ! L - number of encoding levels ! s - interpolation order ! Output: v - encoded cell values ! ! Memory usage in encoding: ! ! Initial data (fine grid values) ! ! ======================================================== ! |v(1,0)|v(2,0)|v(3,0)|v(4,0)|v(5,0)|v(6,0)|v(7,0)|v(8,0) ! ======================================================== ! ! After 1 stage of encoding ! ! ======================================================== ! |v(1,1)|v(2,1)|v(3,1)|v(4,1)|d(1,1)|d(2,1)|d(3,1)|d(4,1) ! ======================================================== ! ! After 2 stages of encoding ! ! ======================================================== ! |v(1,2)|v(2,2)|d(1,2)|d(2,2)|d(1,1)|d(2,1)|d(3,1)|d(4,1) ! ======================================================== ! ! After 3 stages of encoding ! ! ======================================================== ! |v(1,3)|d(1,3)|d(1,2)|d(2,2)|d(1,1)|d(2,1)|d(3,1)|d(4,1) ! ======================================================== ! REAL, ALLOCATABLE, DIMENSION (:) :: vsave, vPrim ALLOCATE(vsave(N)) ALLOCATE(vPrim(0:N)) IF (MOD(N,2**L) .NE. 0) THEN PRINT 1001,N,2**L,L 1001 FORMAT('Error: Number of cells ',i4,' must be a multiple of ',i4, & ' to permit encoding to level ',i4) STOP END IF Nk=N DO k=1,L vsave(1:Nk)=v(1:Nk) Nk=Nk/2 ! Compute average of cell values and values of primitive function at gridpoints ! Note: other averaging procedures might be used here ! such as averaging with a weighting function vPrim(0)=0.0 DO j=1,Nk j2=j+j v(j)=0.5*(vsave(j2-1)+vsave(j2)) vPrim(j)=vPrim(j-1)+v(j) ! Note: hk factor has been left out END DO ! Compute errors in interpolating cell averages from coarser grid DO j=1,Nk j2=j+j v(j+Nk)=vsave(j2-1)-2*(PointInterp(vPrim,j,Nk)-vPrim(j-1)) END DO END DO DEALLOCATE(vsave) DEALLOCATE(vPrim) END SUBROUTINE CellEncode ! Function to decode function cell average values SUBROUTINE CellDecode(N,v,L) INTEGER N,L REAL v(N) ! Input: N - number of cells, must be a multiple of 2**L ! v - encoded cell average values and interpolation errors ! L - number of encoding levels ! s - interpolation order ! Output: v - decoded cell values ! ! Note: On input the coarsest grid cell values and interpolation ! errors are stored in v as specified in the comments on ! memory usage in CellEncode REAL, ALLOCATABLE, DIMENSION (:) :: vsave,vPrim ALLOCATE(vsave(N)) ALLOCATE(vPrim(0:N)) L2=2**L IF (MOD(N,L2) .NE. 0) THEN PRINT 1001,N,L2,L 1001 FORMAT('Error: Number of cells ',i4,' must be a multiple of ',i4, & ' to permit decoding from level ',i4) STOP END IF Nk=N/L2 DO k=L,1,-1 Nk2=Nk+Nk vsave(1:Nk2)=v(1:Nk2) vPrim(0)=0.0 DO j=1,Nk vPrim(j)=vPrim(j-1)+vsave(j) ! Note: hk factor has been left out END DO DO j=1,Nk j2=j+j v(j2-1)=2*(PointInterp(vPrim,j,Nk)-vPrim(j-1))+vsave(j+Nk) v(j2) =2*vsave(j)-v(j2-1) END DO Nk=Nk2 END DO DEALLOCATE(vsave) DEALLOCATE(vPrim) END SUBROUTINE CellDecode ! Function to encode function point values SUBROUTINE PointEncode(N,v,L) INTEGER N,L REAL v(0:N) ! Input: N - number of cells, must be a multiple of 2**L ! v - cell average values of function to be encoded ! L - number of encoding levels ! s - interpolation order ! Output: v - encode cell values ! ! Memory usage in encoding: ! ! Initial data (fine grid values) ! ! =============================================================== ! |v(0,0)|v(1,0)|v(2,0)|v(3,0)|v(4,0)|v(5,0)|v(6,0)|v(7,0)|v(8,0) ! =============================================================== ! ! After 1 stage of encoding ! ! =============================================================== ! |v(0,1)|v(1,1)|v(2,1)|v(3,1)|v(4,1)|d(1,1)|d(2,1)|d(3,1)|d(4,1) ! =============================================================== ! ! After 2 stages of encoding ! ! =============================================================== ! |v(0,2)|v(1,2)|v(2,2)|d(1,2)|d(2,2)|d(1,1)|d(2,1)|d(3,1)|d(4,1) ! =============================================================== ! ! After 3 stages of encoding ! ! =============================================================== ! |v(0,3)|v(1,3)|d(1,3)|d(1,2)|d(2,2)|d(1,1)|d(2,1)|d(3,1)|d(4,1) ! =============================================================== ! REAL, ALLOCATABLE, DIMENSION (:) :: vsave ALLOCATE(vsave(0:N)) IF (MOD(N,2**L) .NE. 0) THEN PRINT 1001,N,2**L,L 1001 FORMAT('Error: Number of points ',i4,' must be a multiple of ',i4, & ' to permit encoding to level ',i4) STOP END IF Nk=N DO k=1,L DO j=0,Nk vsave(j)=v(j) END DO Nk=Nk/2 ! Store even-numbered point values DO j=0,Nk j2=j+j v(j)=vsave(j2) END DO ! Compute errors in interpolating mid-point values from coarser grid DO j=1,Nk j2=j+j v(j+Nk)=vsave(j2-1)-PointInterp(v,j,Nk) END DO END DO DEALLOCATE(vsave) END SUBROUTINE PointEncode ! Function to decode function point values SUBROUTINE PointDecode(N,v,L) INTEGER N,L REAL v(0:N) ! Input: N - number of cells, must be a multiple of 2**L ! v - encoded cell average values and interpolation errors ! L - number of encoding levels ! s - interpolation order ! Output: v - decoded cell values ! ! Note: On input the coarsest grid cell values and interpolation ! errors are stored in v as specified in the comments on ! memory usage in PointEncode REAL, ALLOCATABLE, DIMENSION (:) :: vsave ALLOCATE(vsave(0:N)) L2=2**L IF (MOD(N,L2) .NE. 0) THEN PRINT 1001,N,L2,L 1001 FORMAT('Error: Number of points ',i4,' must be a multiple of ',i4, & ' to permit decoding from level ',i4) STOP END IF DO k=L,1,-1 k2=2**k Nk=N/k2 Nk2=Nk+Nk vsave(0:Nk2)=v(0:Nk2) DO j=1,Nk ! since v(0)=vsave(0) automatically j2=j+j v(j2)=vsave(j) v(j2-1)=PointInterp(vsave,j,Nk)+vsave(j+Nk) END DO Nk=Nk+Nk END DO DEALLOCATE(vsave) END SUBROUTINE PointDecode ! The specific numerical flux computation procedure used is encapsulated ! in this function SUBROUTINE NumericFlux(nFlux,Kflux,mcomp,N,nbc,vn,jstride,fn,j1,j2,sigma) REAL vn(1-nbc:N+nbc,mcomp),fn(0:N,mcomp) REAL f(mcomp,(-Kflux+1):Kflux),v(mcomp,(-Kflux+1):Kflux) REAL sigma jv=j1*jstride ! Numerical fluxes are computed at indices j1,j1+1,...,j2 in the packed fn array ! Finest grid values are always used to compute the fluxes. This is implemented ! by using a stride in the v array which is dependent on the grid level DO j=j1,j2 DO i=-Kflux+1,Kflux ! Compute physical fluxes required by stencil width v(:,i)=vn(jv+i,:) CALL PhysicalFlux(mcomp,v(1,i),f(1,i)) END DO SELECT CASE (nFlux) CASE (1) ! Lax-Friedrichs method fn(j,:)=0.5*((f(:,0)+f(:,1))-(v(:,1)-v(:,0))/sigma) CASE DEFAULT PRINT 1001,nFlux 1001 FORMAT('Numeric flux computation procedure nFlux=',i2,' is not coded yet.') STOP END SELECT jv=jv+jstride END DO fn(0,:)=fn(1,:) END SUBROUTINE NumericFlux ! The physical fluxes from the conservation law SUBROUTINE PhysicalFlux(mcomp,v,f) PARAMETER (gamma=1.4,gamma1=gamma-1.0) REAL v(mcomp),f(mcomp) rho1=1.0/v(1) u=v(2)*rho1 rhou2=v(2)*u p=gamma1*(v(3)-0.5*rhou2) f(1)=v(2) f(2)=p+rhou2 f(3)=u*(p+v(3)) END SUBROUTINE PhysicalFlux ! Carry out explicit time steps of Harten's multiresolution algorithm SUBROUTINE TimeStep(N,nbc,mcomp,vn,delt,delx,tfinal,CFL,nFlux,Kflux,nPlot,p,L,eps) REAL vn(1-nbc:N+nbc,mcomp),eps(L,mcomp),delt,delx,tfinal INTEGER nPlot,p,L LOGICAL, ALLOCATABLE, DIMENSION(:) :: iflag REAL, ALLOCATABLE, DIMENSION (:) :: finterp REAL, ALLOCATABLE, DIMENSION (:,:) :: vnf,vnp1,fbar REAL sigma ! Check that the grids embed properly L2=2**L IF (MOD(N,L2) .NE. 0) THEN PRINT 1001,N,L2,L 1001 FORMAT('Error: Number of points ',i4,' must be a multiple of ',i4, & ' to permit time stepping on ',i4, ' levels') STOP END IF ! Allocate space for flags on whether to use direct evaluation or interpolation ! and work space: vnp1 - new cell averages; fbar - fluxes evaluated directly ! finterp - fluxes evaluated by interpolation ALLOCATE(iFlag(N), vnf(1-nbc:N+nbc,mcomp),vnp1(1-nbc:N+nbc,mcomp), & fbar(0:N,mcomp),finterp(mcomp),STAT=ierror) IF (ierror .NE. 0) THEN PRINT 1002,ierror 1002 FORMAT('Error ',i6,' in allocating iFlag, vnp1, fbar, finterp arrays.') STOP END IF ! Start of time stepping algorithm t=0.0 nStep=0 delt=delt*CFL ! A quick hack. The correct way to go about it is to check maximum ! wave speeds. ! Encode initial approximation DO m=1,mcomp CALL CellEncode(N,vn(1,m),L) END DO DO WHILE (t 2**(p+1)*eps(k,m)) .AND. (k>1)) THEN Nvm1=Nv+Nv joff=j-Nv j2=joff+joff iFlag(Nvm1+j2-1)=.TRUE. iFlag(Nvm1+j2) =.TRUE. END IF END IF END DO END DO Nd=Nd/2 ! go to next coarser level Nv=Nv/2 END DO ! Step (ii) Decode truncated cell values - these are needed for flux computations vnf=vn DO m=1,mcomp CALL CellDecode(N,vnf(1,m),L) END DO ! Step (iii) Evaluate coarsest grid fluxes using the truncated, ! finest grid values L2=2**L NL=N/L2 sigma=delt/delx fbar=0.0 CALL NumericFlux(nFlux,Kflux,mcomp,N,nbc,vnf,L2,fbar,1,NL,sigma) ! Advance v at coarse grid locations sigma=sigma/L2 DO m=1,mcomp DO j=1,NL vnp1(j,m)=vn(j,m)-sigma*(fbar(j,m)-fbar(j-1,m)) END DO END DO ! Step (iv) Compute the representation errors in v for the next time step Nk=N/L2 ! Index of last level-L v element in encoded vn array Nk2=Nk/2 jstride=L2 DO k=L,1,-1 sigma=2*sigma jstride=jstride/2 ! Carry out 1 level of decoding. Interpolated values will be available at ! all grid points of the next finer (k-1) grid DO m=1,mcomp CALL PointDecode(N,fbar(0,m),1) END DO DO j=1,Nk j2=j+j ! Check to see if interpolated fluxes are accurate enough IF (iFlag(Nk+j)) THEN ! No. Must evaluate fluxes directly finterp(:)=fbar(j2-1,:) CALL NumericFlux(nFlux,Kflux,mcomp,N,nbc,vnf,jstride,fbar,j2-1,j2-1,delt/delx) ! Error in interpolating f carries over as an error in v vnp1(Nk+j,:)=vn(Nk+j,:)-sigma*(fbar(j2-1,:)-finterp(:)) ELSE ! Yes. No correction is required. vnp1(Nk+j,:)=0.0 END IF END DO Nk2=Nk Nk=Nk+Nk END DO t=t+delt nStep=nStep+1 ! Impose boundary conditions. Boundary conditions are imposed on the finest grid level. ! This is just a quick hack right now. ! The problem is to be run only for tfinal