C -*- Mode: Fortran; -*- c----------------------------------------------------------------- c Ravi Samtaney c Copyright 2004 c Princeton Plasma Physics Laboratory c All Rights Reserved c----------------------------------------------------------------- c $Log$ c c----------------------------------------------------------------- c For semi-discrete version we need slope fitting in c each cell. This is done using primitive variables c Indices: c vl, vr : 1, nxlocal+1 interior c vl, vr: 2, nxlocal+1 left boundary only c vl, vr: 1, nxlocal right boundary only c vl, vr: 2, nxlocal both boundaries present c dvdx: 0, nxlocal+1 interior c dvdx: 1, nxlocal+1 left boundary c dvdx: 0, nxlocal right boundary c dvdx: 1, nxlocal both boundaries subroutine ConstructLRStatesXi(v,vl,vr) use mesh use properties c implicit none integer:: i,j,k,l,ilo,ihi,direction double precision,intent(in):: v(ixlo:ixhi,nvar) double precision:: vl(1:nxlocal+1,nvar) double precision:: vr(1:nxlocal+1,nvar) double precision:: dvdx(0:nxlocal+1,nvar) c ilo=1+hasxlo ihi=nxlocal+1-hasxhi c c write(6,*) 'CONSTRUCTLRSTATES XI',hasXlo,hasXhi do i=2,nxlocal-1,1 do l=1,nvar,1 dvdx(i,l)=dmin1(dabs(v(i+1,l)-v(i,l))/0.5D0, & dabs(v(i,l)-v(i-1,l))/0.5D0, & dabs(v(i+1,l)-v(i-1,l))/2.0D0)* & dsign(1.D0,v(i+1,l)-v(i-1,l))* & dmax1(dsign(1.D0,(v(i+1,l)-v(i,l)) & *(v(i,l)-v(i-1,l))),0.0D0) dvdx(i,l)=0.D0 enddo enddo c c Left boundary if(hasxlo.eq.1) then i=1 do l=1,nvar,1 dvdx(i,l)=0.D0 c dvdx(i,l)=(v(i+2,l)-v(i,l)) c dvdx(i,l)=dmin1(dabs(v(i+1,l)-v(i,l))/0.5D0, c & dabs(v(i,l)-v(i-1,l))/0.5D0, c & dabs(v(i+1,l)-v(i-1,l))/2.0D0)* c & dsign(1.D0,v(i+1,l)-v(i-1,l))* c & dmax1(dsign(1.D0,(v(i+1,l)-v(i,l)) c & *(v(i,l)-v(i-1,l))),0.0D0) enddo else do i=0,1,1 do l=1,nvar,1 dvdx(i,l)=dmin1(dabs(v(i+1,l)-v(i,l))/0.5D0, & dabs(v(i,l)-v(i-1,l))/0.5D0, & dabs(v(i+1,l)-v(i-1,l))/2.0D0)* & dsign(1.D0,v(i+1,l)-v(i-1,l))* & dmax1(dsign(1.D0,(v(i+1,l)-v(i,l)) & *(v(i,l)-v(i-1,l))),0.0D0) enddo enddo endif c Right boundary if(hasxhi.eq.1) then i=nxlocal do l=1,nvar,1 dvdx(i,l)=0.D0 dvdx(i,l)=(v(i,l)-v(i-1,l))/2.D0 c$$$ dvdx(i,l)=dmin1(dabs(v(i+1,l)-v(i,l))/0.5D0, c$$$ & dabs(v(i,l)-v(i-1,l))/0.5D0, c$$$ & dabs(v(i+1,l)-v(i-1,l))/2.0D0)* c$$$ & dsign(1.D0,v(i+1,l)-v(i-1,l))* c$$$ & dmax1(dsign(1.D0,(v(i+1,l)-v(i,l)) c$$$ & *(v(i,l)-v(i-1,l))),0.0D0) enddo else do i=nxlocal,nxlocal+1,1 do l=1,nvar,1 dvdx(i,l)=dmin1(dabs(v(i+1,l)-v(i,l))/0.5D0, & dabs(v(i,l)-v(i-1,l))/0.5D0, & dabs(v(i+1,l)-v(i-1,l))/2.0D0)* & dsign(1.D0,v(i+1,l)-v(i-1,l))* & dmax1(dsign(1.D0,(v(i+1,l)-v(i,l)) & *(v(i,l)-v(i-1,l))),0.0D0) enddo enddo endif c c dvdx=0.D0 do i=ilo,ihi,1 vl(i,:)=v(i-1,:)+0.5D0*dvdx(i-1,:) vr(i,:)=v(i,:)-0.5D0*dvdx(i,:) enddo c return end subroutine ConstructLRStatesXi c----------------------------------------------------------------- c For semi-discrete version we need slope fitting in c each cell. This is done using primitive variables c Indices: c vl, vr : 1, nylocal+1 interior c vl, vr: 2, nylocal+1 left boundary only c vl, vr: 1, nylocal right boundary only c vl, vr: 2, nylocal both boundaries present c dvdx: 0, nylocal+1 interior c dvdx: 1, nylocal+1 left boundary c dvdx: 0, nylocal right boundary c dvdx: 1, nylocal both boundaries c For many shapes - eta will be periodic subroutine ConstructLRStatesEta(v,vl,vr) use mesh use properties c implicit none integer:: i,j,k,l,ilo,ihi,direction double precision,intent(in):: v(iylo:iyhi,nvar) double precision:: vl(1:nylocal+1,nvar) double precision:: vr(1:nylocal+1,nvar) double precision:: dvdx(0:nylocal+1,nvar) c ilo=1+hasylo ihi=nylocal+1-hasyhi c c write(6,*) 'CONSTRUCTLRSTATES ETA',hasYlo,hasYhi,ilo,ihi c write(6,*) 'CONSTRUCTLRSTATES PERIOD CHECK',v(ilo,8),v(ihi,8) c write(6,*) 'CONSTRUCTLRSTATES PERIOD CHECK',v(ilo-1,8),v(ihi-1,8) c write(6,*) 'CONSTRUCTLRSTATES PERIOD CHECK',v(ilo+1,8),v(ihi+1,8) do i=ilo,ihi-1,1 do l=1,nvar,1 dvdx(i,l)=dmin1(dabs(v(i+1,l)-v(i,l))/0.5D0, & dabs(v(i,l)-v(i-1,l))/0.5D0, & dabs(v(i+1,l)-v(i-1,l))/2.0D0)* & dsign(1.D0,v(i+1,l)-v(i-1,l))* & dmax1(dsign(1.D0,(v(i+1,l)-v(i,l)) & *(v(i,l)-v(i-1,l))),0.0D0) dvdx(i,l)=0.D0 enddo enddo c c Left boundary if(hasylo.eq.1) then i=1 do l=1,nvar,1 dvdx(i,l)=0.D0 c dvdx(i,l)=(v(i+1,l)-v(i,l))/0.5D0 enddo else do i=0,1,1 do l=1,nvar,1 dvdx(i,l)=0.D0 c$$$ dvdx(i,l)=dmin1(dabs(v(i+1,l)-v(i,l))/0.5D0, c$$$ & dabs(v(i,l)-v(i-1,l))/0.5D0, c$$$ & dabs(v(i+1,l)-v(i-1,l))/2.0D0)* c$$$ & dsign(1.D0,v(i+1,l)-v(i-1,l))* c$$$ & dmax1(dsign(1.D0,(v(i+1,l)-v(i,l)) c$$$ & *(v(i,l)-v(i-1,l))),0.0D0) enddo enddo endif c Right boundary if(hasyhi.eq.1) then i=nylocal do l=1,nvar,1 dvdx(i,l)=0.D0 c dvdx(i,l)=(v(i,l)-v(i-1,l))/0.5D0 c$$$ dvdx(i,l)=dmin1(dabs(v(i+1,l)-v(i,l))/0.5D0, c$$$ & dabs(v(i,l)-v(i-1,l))/0.5D0, c$$$ & dabs(v(i+1,l)-v(i-1,l))/2.0D0)* c$$$ & dsign(1.D0,v(i+1,l)-v(i-1,l))* c$$$ & dmax1(dsign(1.D0,(v(i+1,l)-v(i,l)) c$$$ & *(v(i,l)-v(i-1,l))),0.0D0) enddo else do i=nylocal,nylocal+1,1 do l=1,nvar,1 dvdx(i,l)=0.D0 c$$$ dvdx(i,l)=dmin1(dabs(v(i+1,l)-v(i,l))/0.5D0, c$$$ & dabs(v(i,l)-v(i-1,l))/0.5D0, c$$$ & dabs(v(i+1,l)-v(i-1,l))/2.0D0)* c$$$ & dsign(1.D0,v(i+1,l)-v(i-1,l))* c$$$ & dmax1(dsign(1.D0,(v(i+1,l)-v(i,l)) c$$$ & *(v(i,l)-v(i-1,l))),0.0D0) enddo enddo endif c c dvdx=0.D0 do i=ilo,ihi,1 vl(i,:)=v(i-1,:)+0.5D0*dvdx(i-1,:) vr(i,:)=v(i,:)-0.5D0*dvdx(i,:) enddo c return end subroutine ConstructLRStatesEta c-----------------------------------------------------------------