C -*- Mode: Fortran; -*- c----------------------------------------------------------------------- c Ravi Samtaney c Copyright 2004 c Princeton Plasma Physics Laboratory c All Rights Reserved c c----------------------------------------------------------------------- c $Log: inviscidfluxl.F,v $ c Revision 1.3 2005/10/11 22:04:13 samtaney c Cleaned up slightly. c c Revision 1.2 2005/10/02 19:09:32 samtaney c Put in several options. c c Revision 1.1.1.1 2005/09/12 17:28:10 samtaney c Original source. c c Revision 1.3 2005/08/12 20:00:11 samtaney c Added implicit none. c c Revision 1.2 2005/08/11 20:02:11 samtaney c Added commented code to just upwind all except B. c c Revision 1.1 2004/12/04 15:24:44 samtaney c Original source for LF Fluxes. The toroidal component g(psi0)/R is c subtracted. c c----------------------------------------------------------------------- c Computes fluxes using Lax Friedrichs method. c----------------------------------------------------------------------- c Module RusanovFlux c contains subroutine InviscidFluxXi(finv,vx) use mesh use mesh_common use properties use MagAxisToroidalField c implicit none c This subroutine determines the vector of inviscid fluxes. double precision:: finv(1:nxlocal+1,1:nylocal,1:nzlocal,nvar) double precision,intent(in):: & vx(IXLO:IXHI,IYLO:IYHI,IZLO:IZHI,nvar) c double precision,intent(in):: c & beq(IXLO:IXHI,IYLO:IYHI,IZLO:IZHI,3) integer:: direction,ilo,ihi integer:: i,j,k,l c double precision:: vd(IXLO:IXHI,nvar) double precision:: vl(1:nxlocal+1,nvar) double precision:: vr(1:nxlocal+1,nvar) double precision:: wl(1:nxlocal+1,nvar) double precision:: wr(1:nxlocal+1,nvar) double precision:: btt(1:nxlocal+1) double precision:: fn(1:nxlocal+1,nvar) double precision:: fxi(1:nxlocal+1,nvar) double precision:: fc(1:nxlocal+1,nvar) double precision:: vrp(1:nxlocal+1,nvar) double precision:: Rclocal(0:nxlocal+1) double precision:: costheta(1:nxlocal+1) double precision:: sintheta(1:nxlocal+1) double precision:: costhetaC(0:nxlocal+1) double precision:: sinthetaC(0:nxlocal+1) double precision:: dlen(1:nxlocal+1) c ilo=1+hasxlo ihi=nxlocal+1-hasxhi c write(6,*) 'INVISCIDFLUX XI',ilo,ihi, hasxlo do k=1,nzlocal,1 do j=1,nylocal,1 vd(ixlo:ixhi,:)=vx(ixlo:ixhi,j,k,:) btt(ilo:ihi)=Btaxis/RXi(ilo:ihi,j) costheta(ilo:ihi)=dZdEta(ilo:ihi,j) & /dlXi(ilo:ihi,j) sintheta(ilo:ihi)=-dRdEta(ilo:ihi,j) & /dlXi(ilo:ihi,j) call ConstructLRStatesXi(vd,vl,vr) call RemapVariablesToLocalXi(vl,wl,costheta, sintheta, & ilo,ihi) call RemapVariablesToLocalXi(vr,wr,costheta, sintheta, & ilo,ihi) call InviscidFluxRP(fn,vrp,wl,wr,btt,1,nxlocal+1, & ilo,ihi,nvar,1) c call RemapFluxXi(fn,fxi,costheta,sintheta,ilo,ihi) do l=1,nvar,1 finv(ilo:ihi,j,k,l)=RXi(ilo:ihi,j)*fxi(ilo:ihi,l) & *dlXi(ilo:ihi,j) enddo l=7 finv(ilo:ihi,j,k,l)=fxi(ilo:ihi,l) & *dlXi(ilo:ihi,j) c enddo enddo c return end subroutine InviscidFluxXi c c----------------------------------------------------------------------- subroutine RemapVariablesToLocalXi(v,w,costheta,sintheta, & ilo,ihi) use mesh_parms implicit none integer:: ilo,ihi double precision:: v(1:nxlocal+1,nvar) double precision:: w(1:nxlocal+1,nvar) double precision:: costheta(1:nxlocal+1) double precision:: sintheta(1:nxlocal+1) c w(ilo:ihi,1)=v(ilo:ihi,1) w(ilo:ihi,2)=(v(ilo:ihi,2)*costheta(ilo:ihi)+ & v(ilo:ihi,3)*sintheta(ilo:ihi)) w(ilo:ihi,3)=(-v(ilo:ihi,2)*sintheta(ilo:ihi)+ & v(ilo:ihi,3)*costheta(ilo:ihi)) w(ilo:ihi,4)=v(ilo:ihi,4) c w(ilo:ihi,4)=-v(ilo:ihi,4) w(ilo:ihi,5)=(v(ilo:ihi,5)*costheta(ilo:ihi)+ & v(ilo:ihi,6)*sintheta(ilo:ihi)) c KLUDGE c w(ilo:ihi,5)=0.D0 w(ilo:ihi,6)=(-v(ilo:ihi,5)*sintheta(ilo:ihi)+ & v(ilo:ihi,6)*costheta(ilo:ihi)) w(ilo:ihi,7:8)=v(ilo:ihi,7:8) return end subroutine RemapVariablesToLocalXi c----------------------------------------------------------------------- subroutine RemapFluxXi(fn,fxi,costheta,sintheta, & ilo,ihi) use mesh_parms implicit none integer:: ilo,ihi double precision:: fn(1:nxlocal+1,nvar) double precision:: fxi(1:nxlocal+1,nvar) double precision:: costheta(1:nxlocal+1) double precision:: sintheta(1:nxlocal+1) c fxi(ilo:ihi,1)=fn(ilo:ihi,1) fxi(ilo:ihi,2)=fn(ilo:ihi,2)*costheta(ilo:ihi)- & fn(ilo:ihi,3)*sintheta(ilo:ihi) fxi(ilo:ihi,3)=fn(ilo:ihi,2)*sintheta(ilo:ihi)+ & fn(ilo:ihi,3)*costheta(ilo:ihi) fxi(ilo:ihi,4)=fn(ilo:ihi,4) fxi(ilo:ihi,5)=fn(ilo:ihi,5)*costheta(ilo:ihi)- & fn(ilo:ihi,6)*sintheta(ilo:ihi) fxi(ilo:ihi,6)=fn(ilo:ihi,5)*sintheta(ilo:ihi)+ & fn(ilo:ihi,6)*costheta(ilo:ihi) fxi(ilo:ihi,7:8)=fn(ilo:ihi,7:8) c return end subroutine RemapFluxXi c----------------------------------------------------------------------- subroutine InviscidFluxEta(finv,vx) use mesh use mesh_common use properties use MagAxisToroidalField implicit none c c This subroutine determines the vector of inviscid fluxes. double precision:: finv(1:nxlocal,1:nylocal+1,1:nzlocal,nvar) double precision,intent(in):: & vx(IXLO:IXHI,IYLO:IYHI,IZLO:IZHI,nvar) integer:: direction,ilo,ihi integer:: i,j,k,l c double precision:: vd(IYLO:IYHI,nvar) double precision:: vl(1:nylocal+1,nvar) double precision:: vr(1:nylocal+1,nvar) double precision:: wl(1:nylocal+1,nvar) double precision:: wr(1:nylocal+1,nvar) double precision:: fn(1:nylocal+1,nvar) double precision:: feta(1:nylocal+1,nvar) double precision:: fc(1:nylocal+1,nvar) double precision:: vrp(1:nylocal+1,nvar) double precision:: Rclocal(0:nylocal+1) double precision:: costheta(1:nylocal+1) double precision:: sintheta(1:nylocal+1) double precision:: costhetaC(0:nylocal+1) double precision:: sinthetaC(0:nylocal+1) double precision:: btt(1:nylocal+1) c ilo=1+hasylo ihi=nylocal+1-hasyhi c write(6,*) 'INVISCIDFLUX ETA',ilo,ihi do k=1,nzlocal,1 do i=1,nxlocal,1 vd(iylo:iyhi,:)=vx(i,iylo:iyhi,k,:) btt(ilo:ihi)=Btaxis/REta(i,ilo:ihi) costheta(ilo:ihi)=-dZdXi(i,ilo:ihi) & /dlEta(i,ilo:ihi) sintheta(ilo:ihi)=dRdXi(i,ilo:ihi) & /dlEta(i,ilo:ihi) call ConstructLRStatesEta(vd,vl,vr) call RemapVariablesToLocalEta(vl,wl,costheta, sintheta, & ilo,ihi) call RemapVariablesToLocalEta(vr,wr,costheta, sintheta, & ilo,ihi) call InviscidFluxRP(fn,vrp,wl,wr,btt,1,nylocal+1, & ilo,ihi,nvar,2) c call RemapFluxEta(fn,feta,costheta,sintheta,ilo,ihi) do l=1,nvar,1 finv(i,ilo:ihi,k,l)=REta(i,ilo:ihi)*feta(ilo:ihi,l) & *dlEta(i,ilo:ihi) enddo l=7 finv(i,ilo:ihi,k,l)=feta(ilo:ihi,l) & *dlEta(i,ilo:ihi) enddo enddo c return end subroutine InviscidFluxEta c c----------------------------------------------------------------------- subroutine RemapVariablesToLocalEta(v,w,costheta,sintheta & ,ilo,ihi) use mesh_parms implicit none integer:: ilo,ihi double precision:: v(1:nylocal+1,nvar) double precision:: w(1:nylocal+1,nvar) double precision:: costheta(1:nylocal+1) double precision:: sintheta(1:nylocal+1) c w(ilo:ihi,1)=v(ilo:ihi,1) w(ilo:ihi,2)=v(ilo:ihi,2)*costheta(ilo:ihi)+ & v(ilo:ihi,3)*sintheta(ilo:ihi) w(ilo:ihi,3)=-v(ilo:ihi,2)*sintheta(ilo:ihi)+ & v(ilo:ihi,3)*costheta(ilo:ihi) w(ilo:ihi,4)=v(ilo:ihi,4) c w(ilo:ihi,4)=-v(ilo:ihi,4) w(ilo:ihi,5)=v(ilo:ihi,5)*costheta(ilo:ihi)+ & v(ilo:ihi,6)*sintheta(ilo:ihi) w(ilo:ihi,6)=-v(ilo:ihi,5)*sintheta(ilo:ihi)+ & v(ilo:ihi,6)*costheta(ilo:ihi) cKLUDGE c w(ilo:ihi,6)=0.D0 w(ilo:ihi,7)=v(ilo:ihi,7) w(ilo:ihi,8)=v(ilo:ihi,8) return end subroutine RemapVariablesToLocalEta c----------------------------------------------------------------------- subroutine RemapFluxEta(fn,feta,costheta,sintheta,ilo,ihi) use mesh integer:: ilo,ihi double precision:: fn(1:nylocal+1,nvar) double precision:: feta(1:nylocal+1,nvar) double precision:: costheta(1:nylocal+1) double precision:: sintheta(1:nylocal+1) c feta(ilo:ihi,1)=fn(ilo:ihi,1) feta(ilo:ihi,2)=fn(ilo:ihi,2)*costheta(ilo:ihi)- & fn(ilo:ihi,3)*sintheta(ilo:ihi) feta(ilo:ihi,3)=fn(ilo:ihi,2)*sintheta(ilo:ihi)+ & fn(ilo:ihi,3)*costheta(ilo:ihi) feta(ilo:ihi,4)=fn(ilo:ihi,4) feta(ilo:ihi,5)=fn(ilo:ihi,5)*costheta(ilo:ihi)- & fn(ilo:ihi,6)*sintheta(ilo:ihi) feta(ilo:ihi,6)=fn(ilo:ihi,5)*sintheta(ilo:ihi)+ & fn(ilo:ihi,6)*costheta(ilo:ihi) feta(ilo:ihi,7:8)=fn(ilo:ihi,7:8) c feta(ilo:ihi,8)=fn(ilo:ihi,8) return end subroutine RemapFluxEta c----------------------------------------------------------------------- subroutine InviscidFluxPhi(finv,vx) use mesh use mesh_common use properties use MagAxisToroidalField implicit none c c This subroutine determines the vector of inviscid fluxes. double precision:: finv(1:nxlocal,1:nylocal,1:nzlocal+1,nvar) double precision,intent(in):: & vx(IXLO:IXHI,IYLO:IYHI,IZLO:IZHI,nvar) integer:: direction,ilo,ihi integer:: i,j,k,l c double precision:: vd(IZLO:IZHI,nvar) double precision:: vl(1:nzlocal+1,nvar) double precision:: vr(1:nzlocal+1,nvar) double precision:: wl(1:nzlocal+1,nvar) double precision:: wr(1:nzlocal+1,nvar) double precision:: fn(1:nzlocal+1,nvar) double precision:: fphi(1:nzlocal+1,nvar) double precision:: fc(1:nzlocal+1,nvar) double precision:: vrp(1:nzlocal+1,nvar) double precision:: btt(1:nzlocal+1) c ilo=1 ihi=nzlocal+1 c write(6,*) 'INVISCIDFLUX PHI',ilo,ihi do j=1,nylocal,1 do i=1,nxlocal,1 vd(izlo:izhi,:)=vx(i,j,izlo:izhi,:) btt(ilo:ihi)=Btaxis/Rc(i,j) call ConstructLRStatesPhi(vd,vl,vr) call RemapVariablesToLocalPhi(vl,wl,ilo,ihi) call RemapVariablesToLocalPhi(vr,wr,ilo,ihi) call InviscidFluxRP(fn,vrp,wl,wr,btt,1,nzlocal+1, & ilo,ihi,nvar,3) call RemapFluxPhi(fn,fphi,ilo,ihi) finv(i,j,ilo:ihi,:)=fphi(ilo:ihi,:) enddo enddo c return end subroutine InviscidFluxPhi c c----------------------------------------------------------------------- subroutine RemapVariablesToLocalPhi(v,w,ilo,ihi) use mesh implicit none integer:: ind integer:: ilo,ihi double precision:: v(1:nzlocal+1,nvar) double precision:: w(1:nzlocal+1,nvar) c w(ilo:ihi,1)=v(ilo:ihi,1) w(ilo:ihi,2)=v(ilo:ihi,4) c w(ilo:ihi,2)=-v(ilo:ihi,4) w(ilo:ihi,3)=v(ilo:ihi,2) w(ilo:ihi,4)=v(ilo:ihi,3) w(ilo:ihi,5)=v(ilo:ihi,7) w(ilo:ihi,6)=v(ilo:ihi,5) w(ilo:ihi,7)=v(ilo:ihi,6) w(ilo:ihi,8)=v(ilo:ihi,8) return end subroutine RemapVariablesToLocalPhi c----------------------------------------------------------------------- subroutine RemapFluxPhi(fn,fphi,ilo,ihi) use mesh implicit none integer:: ilo,ihi double precision:: fn(1:nzlocal+1,nvar) double precision:: fphi(1:nzlocal+1,nvar) c fphi(ilo:ihi,1)=fn(ilo:ihi,1) fphi(ilo:ihi,4)=fn(ilo:ihi,2) fphi(ilo:ihi,2)=fn(ilo:ihi,3) fphi(ilo:ihi,3)=fn(ilo:ihi,4) fphi(ilo:ihi,7)=fn(ilo:ihi,5) fphi(ilo:ihi,5)=fn(ilo:ihi,6) fphi(ilo:ihi,6)=fn(ilo:ihi,7) fphi(ilo:ihi,8)=fn(ilo:ihi,8) return end subroutine RemapFluxPhi c----------------------------------------------------------------------- subroutine InviscidFluxRP(finv,vm,vl,vr,btt,inlo, & inhi,ilo,ihi,nvar,direction) use properties use MagAxisToroidalField use alphas use evalues implicit none integer:: inlo,inhi, ilo,ihi,nvar double precision:: finv(inlo:inhi,nvar) double precision,intent(in):: vl(inlo:inhi,nvar) double precision,intent(in):: vr(inlo:inhi,nvar) double precision,intent(in):: btt(inlo:inhi) double precision:: vm(inlo:inhi,nvar) double precision:: rho,u,v,w,bx,by,bz double precision:: eint, lam, press double precision:: alamda(ILO:IHI,2) double precision:: durl(ILO:IHI,nvar) double precision:: finvl(nvar) double precision:: finvr(nvar) double precision:: bxsq,bysq,bzsq integer:: i,j,k,l,m,direction,ii double precision:: Btlocal double precision:: lamdaMax c c This subroutine determines the vector of inviscid fluxes. c c v: Vector of variables - rho, u,v,w,Bx,By,Bz,p c rho=density, u=x-velocity, v=y-velocity c w=z-velocity c Bx,By,Bz= Magnetic Field c p=total pressure pressure c c vl(i) = Left state after linear representation in cell, i . c at the i-th interface c vr(i) = Right state after linear representation in cell, i+1 . c at the i-th interface c c double precision:: utilde(ilo:ihi,nvar) c utilde = vector of Roe averaged variables. c utilde(1,2,3,4,5,6,7,8) = {rhotilde, utilde, vtilde, wtilde, c Bxtilde,Bytilde,Bztilde,ptilde} c c Set the Roe averaged variables depending on vl,vr call SetRoeVariables(utilde,vl,vr,inlo,inhi,ilo,ihi) c alamda=0.D0 durl=0.D0 c call SetEigenValues(utilde,alamda,ilo,ihi, & btt,inlo,inhi,direction,nvar) call SetDurl(durl,vl,vr,inlo,inhi,ilo,ihi) c do i=ilo,ihi,1 lamdaMax=maxval(dabs(alamda(i,1:8))) c rho=vl(i,1) u=vl(i,2); v=vl(i,3); w=vl(i,4) bx=vl(i,5); by=vl(i,6); bz=vl(i,7) bxsq=bx*bx bysq=by*by bzsq=bz*bz press=vl(i,8)-0.5D0*(bxsq+bysq+bzsq) finvl(1)=rho*u finvl(2)=rho*u*u +press+0.5D0*(by*by+bz*bz-bx*bx) finvl(3)=rho*u*v-bx*by finvl(4)=rho*u*w-bx*bz eint=press/(gamma-1.D0) finvl(5)=0.D0 finvl(6)=u*by-v*bx finvl(7)=u*(bz)-w*bx finvl(8)=(0.5D0*rho*(u*u+v*v+w*w)+eint+ & press+(bx*bx+by*by+bz*bz))*vl(i,2)- & bx*(u*bx+v*by+w*bz) c if(direction.eq.1.or.direction.eq.2) then finvl(2)=finvl(2)+Btt(i)*bz finvl(7)=finvl(7)+Btt(i)*u finvl(8)=finvl(8)+Btt(i)*bz*u else if(direction.eq.3) then finvl(3)=finvl(3)-Btt(i)*by finvl(4)=finvl(4)-Btt(i)*bz finvl(6)=finvl(6)-Btt(i)*v finvl(7)=finvl(7)-Btt(i)*w finvl(8)=finvl(8)-Btt(i)*(by*v+bz*w) endif c rho=vr(i,1) u=vr(i,2); v=vr(i,3); w=vr(i,4) bx=vr(i,5); by=vr(i,6); bz=vr(i,7) bxsq=bx*bx bysq=by*by bzsq=bz*bz press=vr(i,8)-0.5D0*(bxsq+bysq+bzsq) finvr(1)=rho*u finvr(2)=rho*u*u +press+0.5D0*(by*by+bz*bz-bx*bx) finvr(3)=rho*u*v-bx*by finvr(4)=rho*u*w-bx*bz eint=press/(gamma-1.D0) finvr(5)=0.D0 finvr(6)=u*by-v*bx finvr(7)=u*(bz)-w*bx finvr(8)=(0.5D0*rho*(u*u+v*v+w*w)+eint+ & press+(bx*bx+by*by+bz*bz))*vr(i,2)- & bx*(u*bx+v*by+w*bz) c c if(direction.eq.1.or.direction.eq.2) then finvr(2)=finvr(2)+Btt(i)*bz finvr(7)=finvr(7)+Btt(i)*u finvr(8)=finvr(8)+Btt(i)*bz*u else if(direction.eq.3) then finvr(3)=finvr(3)-Btt(i)*by finvr(4)=finvr(4)-Btt(i)*bz finvr(6)=finvr(6)-Btt(i)*v finvr(7)=finvr(7)-Btt(i)*w finvr(8)=finvr(8)-Btt(i)*(by*v+bz*w) endif c Do not upwind B c$$$ finv(i,1:4)=0.5D0*(finvl(1:4)+finvr(1:4)) c$$$ & -0.5D0*lamdaMax*durl(i,1:4) c$$$ finv(i,5)=0.D0 c$$$ finv(i,6:7)=0.5D0*(finvl(6:7)+finvr(6:7)) c$$$ finv(i,8)=0.5D0*(finvl(8)+finvr(8)) c$$$ & -0.5D0*lamdaMax*durl(i,8) c Uncomment lines below for upwinding all finv(i,1:8)=0.5D0*(finvl(1:8)+finvr(1:8)) & -0.5D0*lamdaMax*durl(i,:) finv(i,5)=0.D0 enddo c c return end subroutine InviscidFluxRP c end Module RusanovFlux