!======================================================================== ! BEARCLAW Boundary Embedded Adaptive Refinement Conservation LAW package !======================================================================== ! (c) Copyright Sorin Mitran, 2000 ! Department of Applied Mathematics ! University of Washington ! mitran@amath.washington.edu ! ----------------------------------------------------------------------- ! This code may be freely used for educational and research purposes. ! For any other use please contact the author. ! ----------------------------------------------------------------------- ! File: classicbear.f90 ! Purpose: Coding of some classic schemes as the time stepping ! routine for BEARCLAW ! Contains: ! Revision History: Ver. 1.0 Aug. 2001 Sorin Mitran ! ----------------------------------------------------------------------- ! Comments: ! 1) This is an example implementation of the BEARez TimeStep function. ! The basic task a TimeStep function must perform is to provide an ! updated Info%qnew reflecting field values at t+dt, starting ! from the field values at time t given in Info%qold. ! ----------------------------------------------------------------------- !************ MODULE scheme !************ ! User definitions of data structures associated with each node USE NodeInfoDef ! Problem specific routines (setprob,qinit,b4step,src,physflux) USE problem IMPLICIT NONE SAVE PRIVATE ! Everything is implicitly PRIVATE PUBLIC numflux,updateq ! Global variables INTEGER, PARAMETER :: LaxFriedrichs=1,Upwind=2 INTEGER, PARAMETER :: LaxWendroff=3,BeamWarming=4,Fromm=5,minmod=6,superbee=7,MC=8,vanLeer=9 REAL (KIND=xPrec), PARAMETER :: zero=0.0 !------- CONTAINS !------- ! Numerical flux routine: ! Compute the numerical flux along dimension n ! Slice is at indices indx SUBROUTINE numflux(n,indx,inlow,inhigh,Info) ! Interface declarations INTEGER n,indx(MaxDims),inlow,inhigh TYPE (NodeInfo) :: Info ! Internal declarations INTEGER mbc,mx,m,j,iError,nq,i1,i2,i,indx2(MaxDims),nDim INTEGER, ALLOCATABLE, DIMENSION(:) :: Ip1,Im1 REAL (KIND=qPrec) :: fact,dtdx,dtdx2 REAL (KIND=qPrec), POINTER, DIMENSION (:) :: dtdx1D REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: q1D,fql,fqr,sl,sr,numfql,numfqr,qcenter,qedge REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: fq,s,work REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: Ar,Al,A,Acenter,Aedge ! Associate local names with fields of grid data structure q1D=>Info%q1D; dtdx1D=>Info%dtdx1D fql=>Info%fql; numfql=>Info%numfql; sl=>Info%sl; Al=>Info%Aql fqr=>Info%fqr; numfqr=>Info%numfqr; sr=>Info%sr; Ar=>Info%Aqr ! More suggestive names for algorithms which don't differentiate ! between right and left edges of a cell fq=>Info%fqr; s=>Info%sr; A=>Info%Aqr; ! Clean out memory contents fql=zero; numfql=zero; sl=zero; Al=zero fqr=zero; numfqr=zero; sr=zero; Ar=zero ! Current number of interior cellsm ghost cells, field variables mx=Info%mXnow; mbc=Info%mbc; nq=Info%NrVars ! Algorithms are programmed for uniform spacing dtdx=Info%dt/Info%dX(n); dtdx2=0.5*dtdx ! Useful index range vectors ALLOCATE(Ip1(nq),Im1(nq),STAT=iError) IF (iError /= 0) THEN PRINT *,'Cannot allocate index range arrays in numflux' STOP END IF Im1=SPREAD(-1,1,nq); Ip1=SPREAD(1,1,nq) SELECT CASE (Info%method(8)) CASE (LaxFriedrichs) fact=1.0/2.0**Info%nDim/dtdx CALL extract1D(n,indx,inlow,inhigh,Info,q1D) CALL physflux(n,indx,RequestFluxes,Info,q1D,fq,s,A) numfqr = 0.5*(fq + EOSHIFT(fq,Ip1)) - fact*(EOSHIFT(q1D,Ip1) - q1D) numfql = 0.5*(EOSHIFT(fq,Im1) + fq) - fact*(q1D - EOSHIFT(q1D,Im1)) CASE (Upwind) CALL extract1D(n,indx,inlow,inhigh,Info,q1D) CALL physflux(n,indx,RequestFluxes,Info,q1D,fq,s,A) CALL physflux(n,indx,RequestSpeeds,Info,q1D,fq,s,A) ! Numerical flux at left edge of cell WHERE (szero) numfqr = fq ELSEWHERE numfqr = EOSHIFT(fq,Ip1) END WHERE CASE (LaxWendroff) ! More descriptive names Acenter=>Info%Aql; Aedge=>Info%Aqr; qcenter=>Info%q1D; qedge=>Info%sr ! We won't need the speeds, so use the memory space for ! evaluating cell edge values work=>Info%sl qedge=0.; Aedge=0. qcenter=0.; Acenter=0. work=0. ! Get the physical fluxes at the cell centers CALL extract1D(n,indx,inlow,inhigh,Info,qcenter) CALL physflux(n,indx,RequestFluxes,Info,qcenter,fq,s,Acenter) ! Get the flux Jacobians at the cell centers CALL physflux(n,indx,RequestJacobians,Info,qcenter,fq,s,Acenter) numfqr=0.5*(fq+EOSHIFT(fq,Ip1)) numfql=0.5*(EOSHIFT(fq,Im1)+fq) ! Evaluate cell edge field values by arithmetic average qedge=0.5*(qcenter+EOSHIFT(qcenter,Ip1)) CALL physflux(n,indx,RequestJacobians,Info,qedge,fq,s,Aedge) fact=0.5*dtdx work=EOSHIFT(fq,Ip1)-fq DO i=2-mbc,mx+mbc-1 work(i,1:nq)=fact*MATMUL(Aedge(i,1:nq,1:nq),work(i,1:nq)) END DO numfqr=numfqr-work numfql=numfql-EOSHIFT(work,Im1) ! Include cross-derivative terms DO nDim=1,Info%nDim indx2=indx IF (nDim==n) CYCLE fact=0.125*Info%dt/Info%dX(nDim) indx2(nDim)=indx(nDim)+1 CALL extract1D(n,indx2,inlow,inhigh,Info,qcenter) CALL physflux(n,indx,RequestFluxes,Info,qcenter,fq,s,Acenter) work=fq indx2(nDim)=indx(nDim)-1 CALL extract1D(n,indx2,inlow,inhigh,Info,qcenter) CALL physflux(n,indx,RequestFluxes,Info,qcenter,fq,s,Acenter) work=work-fq DO i=2-mbc,mx+mbc-1 work(i,1:nq)=fact*MATMUL(Acenter(i,1:nq,1:nq),work(i,1:nq)) END DO numfqr=numfqr-work-EOSHIFT(work,Ip1) numfql=numfql-EOSHIFT(work,Im1)-work END DO CASE (BeamWarming) CASE (Fromm) CASE (minmod) CASE (superbee) CASE (MC) CASE (vanLeer) CASE DEFAULT PRINT *,'Bad method(8)=',Info%method(8),'. No such scheme implemented.' STOP END SELECT DEALLOCATE(Im1,Ip1) END SUBROUTINE numflux ! Update field variable routine SUBROUTINE updateq(n,indx,inlow,inhigh,Info) ! Interface declarations INTEGER n,indx(4),inlow,inhigh TYPE (NodeInfo) :: Info ! Internal declarations INTEGER mbc,mcapa,i1,i2,i3,i4,i1Dlow,i1Dhigh,mX,nq REAL (KIND=qPrec), POINTER, DIMENSION (:) :: dtdx1D mcapa=Info%method(6); mbc=Info%mbc; nq=Info%NrVars dtdx1D=>Info%dtdx1D ! Link between standard 1-mbc:mx+mbc index space and AMR index space i1Dlow=inlow+mbc; i1Dhigh=inhigh-mbc i1=indx(1); i2=indx(2); i3=indx(3); i4=indx(4); mX=Info%mXnow SELECT CASE (n) CASE (1) Info%q(i1Dlow:i1Dhigh,i2,i3,i4,1:nq) = Info%q(i1Dlow:i1Dhigh,i2,i3,i4,1:nq) - & SPREAD(dtdx1D(1:mX),1,nq)*(Info%numfqr(1:mX,1:nq)-Info%numfql(1:mX,1:nq)) CASE (2) Info%q(i1,i1Dlow:i1Dhigh,i3,i4,1:nq) = Info%q(i1,i1Dlow:i1Dhigh,i3,i4,1:nq) - & SPREAD(dtdx1D(1:mX),1,nq)*(Info%numfqr(1:mX,1:nq)-Info%numfql(1:mX,1:nq)) CASE (3) Info%q(i1,i2,i1Dlow:i1Dhigh,i4,1:nq) = Info%q(i1,i2,i1Dlow:i1Dhigh,i4,1:nq) - & SPREAD(dtdx1D(1:mX),1,nq)*(Info%numfqr(1:mX,1:nq)-Info%numfql(1:mX,1:nq)) CASE (4) Info%q(i1,i2,i3,i1Dlow:i1Dhigh,1:nq) = Info%q(i1,i2,i3,i1Dlow:i1Dhigh,1:nq) - & SPREAD(dtdx1D(1:mX),1,nq)*(Info%numfqr(1:mX,1:nq)-Info%numfql(1:mX,1:nq)) END SELECT END SUBROUTINE updateq ! Extract a 1D slice of data SUBROUTINE extract1D(n,indx,inlow,inhigh,Info,q1D) ! Interface declarations INTEGER n,indx(4),inlow,inhigh TYPE (NodeInfo) :: Info REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: q1D ! Internal declarations INTEGER i1Dlow,i1Dhigh,i1,i2,i3,i4 ! Transform 1D slice of index space to standard 1-mbc:mx+mbc i1Dlow=1-Info%mbc; i1Dhigh=inhigh-inlow+i1Dlow ! Copy a 1D slice of data along dimension n from qold to q1D i1=indx(1); i2=indx(2); i3=indx(3); i4=indx(4) SELECT CASE (n) CASE (1) q1D(i1Dlow:i1Dhigh,:)=Info%qold(inlow:inhigh,i2,i3,i4,:) CASE (2) q1D(i1Dlow:i1Dhigh,:)=Info%qold(i1,inlow:inhigh,i3,i4,:) CASE (3) q1D(i1Dlow:i1Dhigh,:)=Info%qold(i1,i2,inlow:inhigh,i4,:) CASE (4) q1D(i1Dlow:i1Dhigh,:)=Info%qold(i1,i2,i3,inlow:inhigh,:) END SELECT ! Compute step ratio Info%dtdx1D(i1Dlow:i1Dhigh)=Info%dt/Info%dX(n) END SUBROUTINE extract1D ! End of module CONTAINS clause !**************** END MODULE scheme !****************