* remove support for "General" EOS interface (EOSG_*)

* test suites pass with Intel 11 (no optimizaton) on bethe 
  and gcc 4.4.5




git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@225 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

1 files changed, 0 insertions, 353 deletions

diff --git a/src/GRHydro_RiemannSolveM.F90 b/src/GRHydro_RiemannSolveM.F90
index 8944b26..bd16778 100644
--- a/src/GRHydro_RiemannSolveM.F90
+++ b/src/GRHydro_RiemannSolveM.F90
@@ -140,358 +140,5 @@ end subroutine RiemannSolvePolytypeM
 
 
 
-/*@@
-@routine    RiemannSolveGeneralM
-@date       Tue Mar 19 11:40:20 2002
-@author     Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
-@desc 
-The Riemann solvers for the new general EOS routines.
-This sets the fluxes from the left and right reconstructed
-states, so that after this routine they are effectively 
-scratch space.
-@enddesc 
-@calls     
-@calledby   
-@history 
-
-@endhistory 
-
-@@*/
-
-
-subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS)
-  
-  USE GRHydro_Scalars
-  
-  implicit none
-  
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_PARAMETERS
-  DECLARE_CCTK_FUNCTIONS
-
-  CCTK_INT :: i,j,k, ierr
-  CCTK_REAL, dimension(8) :: tmp_flux, cons_p, cons_m
-  CCTK_REAL, dimension(6) :: prim_p, prim_m
-  CCTK_REAL :: avg_det, avg_alp, avg_beta
-  CCTK_REAL :: gxxh, gyyh, gzzh, gxyh, gxzh, gyzh
-  CCTK_REAL :: avg_betax, avg_betay, avg_betaz
-  CCTK_REAL :: vxtp,vytp,vztp,vxtm,vytm,vztm,ab0p,ab0m,b2p,b2m,bdotvp,bdotvm
-  CCTK_REAL :: wp,wm,v2p,v2m,bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,vA2m,vA2p
-  CCTK_REAL :: Bvecxlowp,Bvecxlowm,Bvecylowp,Bvecylowm,Bveczlowp,Bveczlowm
-  CCTK_REAL :: pressstarp,pressstarm,rhoenth_p,rhoenth_m
-  CCTK_REAL :: velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm
-
-  CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm
-
-  densflux = 0.d0
-  sxflux   = 0.d0
-  syflux   = 0.d0
-  szflux   = 0.d0
-  tauflux  = 0.d0
-  Bvecxflux = 0.d0
-  Bvecyflux = 0.d0
-  Bveczflux = 0.d0
-  if(clean_divergence.ne.0) then
-     psidcflux = 0.d0
-  endif
-
-!!$  Do the EOS call to set the pressure, derivative and cs2
-
-  ierr = EOS_SetGFs(cctkGH, EOS_RiemannCallPlus)
-  ierr = EOS_SetGFs(cctkGH, EOS_RiemannCallMinus)
-
-
-  if (CCTK_EQUALS(riemann_solver,"HLLE")) then
-     
-     call GRHydro_HLLEGeneralM(CCTK_PASS_FTOF)
-     
-     if (evolve_tracer .ne. 0) then
-        
-!!$ No b-field component for tracers!
-        call GRHydro_HLLE_TracerGeneral(CCTK_PASS_FTOF)
-        
-     end if
-     
-  else
-     
-     do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
-        do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
-           do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil        
-              
-!!$        Set the left (p for plus) and right (m for minus, i+1) states
-              
-              cons_p(1)   = densplus(i,j,k) 
-              cons_p(2)   = sxplus(i,j,k)
-              cons_p(3)   = syplus(i,j,k)
-              cons_p(4)   = szplus(i,j,k)
-              cons_p(5)   = tauplus(i,j,k)
-              cons_p(6)   = Bvecxplus(i,j,k)
-              cons_p(7)   = Bvecyplus(i,j,k)
-              cons_p(8)   = Bveczplus(i,j,k)
-              
-              cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset)
-              cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset)
-              cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset)
-              cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset)
-              cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset) 
-              cons_m(6) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset)
-              cons_m(7) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset)
-              cons_m(8) = Bveczminus(i+xoffset,j+yoffset,k+zoffset)
-              
-              prim_p(1)   = rhoplus(i,j,k) 
-              prim_p(2)   = velxplus(i,j,k)
-              prim_p(3)   = velyplus(i,j,k) 
-              prim_p(4)   = velzplus(i,j,k)
-              prim_p(5)   = epsplus(i,j,k)
-              prim_p(6)   = pressplus(i,j,k)
-              rhoenth_p    = prim_p(1)*(1.0d0+prim_p(5))+prim_p(6)
-              
-              prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset)
-              prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset)
-              prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset)
-              prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset)
-              prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) 
-              prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) 
-              rhoenth_m    = prim_m(1)*(1.0d0+prim_m(5))+prim_m(6)
-
-              if(clean_divergence.ne.0) then
-                 psidcp = psidcplus(i,j,k)
-                 psidcm = psidcminus(i+xoffset,j+yoffset,k+zoffset)
-              endif
-             
-!!$        Set metric terms at interface
-              
-              if (shift_state .ne. 0) then
-                 avg_betax = 0.5d0*(betax(i+xoffset,j+yoffset,k+zoffset)+betax(i,j,k))
-                 avg_betay = 0.5d0*(betay(i+xoffset,j+yoffset,k+zoffset)+betay(i,j,k))
-                 avg_betaz = 0.5d0*(betaz(i+xoffset,j+yoffset,k+zoffset)+betaz(i,j,k))
-                 if (flux_direction == 1) then
-                    avg_beta = avg_betax
-                 else if (flux_direction == 2) then
-                    avg_beta = avg_betay
-                 else if (flux_direction == 3) then
-                    avg_beta = avg_betaz
-                 else
-                    call CCTK_WARN(0, "Flux direction not x,y,z")
-                 end if
-              else
-                 avg_beta = 0.d0
-                 avg_betax = 0.d0
-                 avg_betay = 0.d0
-                 avg_betaz = 0.d0
-              end if
-              
-              avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
-              
-              gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + &
-                   gxx(i,j,k))
-              gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + &
-                   gxy(i,j,k))
-              gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + &
-                   gxz(i,j,k))
-              gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + &
-                   gyy(i,j,k))
-              gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + &
-                   gyz(i,j,k))
-              gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + &
-                   gzz(i,j,k))
-              
-              avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)
-              
-              vxtp = prim_p(2)-avg_betax/avg_alp
-              vytp = prim_p(3)-avg_betay/avg_alp
-              vztp = prim_p(4)-avg_betaz/avg_alp
-              vxtm = prim_m(2)-avg_betax/avg_alp
-              vytm = prim_m(3)-avg_betay/avg_alp
-              vztm = prim_m(4)-avg_betaz/avg_alp
-              
-              call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
-                   prim_p(2),prim_p(3),prim_p(4),cons_p(6),cons_p(7),cons_p(8), &
-                   velxlowp,velylowp,velzlowp,Bvecxlowp,Bvecylowp,Bveczlowp, &
-                   bdotvp,b2p,v2p,wp,bxlowp,bylowp,bzlowp)
-              call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
-                   prim_m(2),prim_m(3),prim_m(4),cons_m(6),cons_m(7),cons_m(8), &
-                   velxlowm,velylowm,velzlowm,Bvecxlowm,Bvecylowm,Bveczlowm, &
-                   bdotvm,b2m,v2m,wm,bxlowm,bylowm,bzlowm)
-              
-              ab0p = wp*bdotvp
-              ab0m = wm*bdotvm
-              
-              vA2p = b2p/(rhoenth_p+b2p)
-              vA2m = b2m/(rhoenth_m+b2m)
-              
-!!$ p^*  = p+b^2/2  in Anton et al.
-              pressstarp = prim_p(6)+0.5d0*b2p
-              pressstarm = prim_m(6)+0.5d0*b2m        
-              
-              if (flux_direction == 1) then
-                 
-                 call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),&
-                      cons_p(6),cons_p(7),cons_p(8),&
-                      tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5), &
-                      tmp_flux(6),tmp_flux(7),tmp_flux(8), &
-                      vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, &
-                      avg_det,avg_alp,avg_beta)
-
-                 densflux(i,j,k) = 0.5d0 * tmp_flux(1)
-                 sxflux(i,j,k)   = 0.5d0 * tmp_flux(2)
-                 syflux(i,j,k)   = 0.5d0 * tmp_flux(3)
-                 szflux(i,j,k)   = 0.5d0 * tmp_flux(4)
-                 tauflux(i,j,k)  = 0.5d0 * tmp_flux(5)
-                 Bvecxflux(i,j,k)   = 0.5d0 * tmp_flux(6)
-                 Bvecyflux(i,j,k)   = 0.5d0 * tmp_flux(7)
-                 Bveczflux(i,j,k)   = 0.5d0 * tmp_flux(8)
-
-                 if(clean_divergence.ne.0) then
-                    psidcf = cons_p(6)
-                    psidcflux(i,j,k) = 0.5d0 * psidcf
-                 endif
-                 
-                 call num_x_fluxM(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
-                      cons_m(6),cons_m(7),cons_m(8),&
-                      tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5),&
-                      tmp_flux(6),tmp_flux(7),tmp_flux(8),&
-                      vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, &
-                      avg_det,avg_alp,avg_beta)
-
-                 densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1)
-                 sxflux(i,j,k)   = sxflux(i,j,k)   + 0.5d0 * tmp_flux(2)
-                 syflux(i,j,k)   = syflux(i,j,k)   + 0.5d0 * tmp_flux(3)
-                 szflux(i,j,k)   = szflux(i,j,k)   + 0.5d0 * tmp_flux(4)
-                 tauflux(i,j,k)  = tauflux(i,j,k)  + 0.5d0 * tmp_flux(5)
-                 Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6)
-                 Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
-                 Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
-                 
-                 if(clean_divergence.ne.0) then
-                    psidcf = cons_m(6)
-                    psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
-                 endif
-
-              else if (flux_direction == 2) then
-                 
-                 call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),&
-                      cons_p(7),cons_p(8),cons_p(6),&
-                      tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5),&
-                      tmp_flux(7),tmp_flux(8),tmp_flux(6),&
-                      vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, &
-                      avg_det,avg_alp,avg_beta)
-                 
-                 densflux(i,j,k) = 0.5d0 * tmp_flux(1)
-                 sxflux(i,j,k)   = 0.5d0 * tmp_flux(2)
-                 syflux(i,j,k)   = 0.5d0 * tmp_flux(3)
-                 szflux(i,j,k)   = 0.5d0 * tmp_flux(4)
-                 tauflux(i,j,k)  = 0.5d0 * tmp_flux(5)
-                 Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6)
-                 Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
-                 Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
-                 
-                 if(clean_divergence.ne.0) then
-                    psidcf = cons_p(7)
-                    psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
-                 endif
-
-                 call num_x_fluxM(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
-                      cons_m(7),cons_m(8),cons_m(6),&
-                      tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5),&
-                      tmp_flux(7),tmp_flux(8),tmp_flux(6),&
-                      vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, &
-                      avg_det,avg_alp,avg_beta)
-             
-                 densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1)
-                 sxflux(i,j,k)   = sxflux(i,j,k)   + 0.5d0 * tmp_flux(2)
-                 syflux(i,j,k)   = syflux(i,j,k)   + 0.5d0 * tmp_flux(3)
-                 szflux(i,j,k)   = szflux(i,j,k)   + 0.5d0 * tmp_flux(4)
-                 tauflux(i,j,k)  = tauflux(i,j,k)  + 0.5d0 * tmp_flux(5)
-                 Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6)
-                 Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
-                 Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
-                 
-                 if(clean_divergence.ne.0) then
-                    psidcf = cons_m(7)
-                    psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
-                 endif
-
-              else if (flux_direction == 3) then
-                 
-                 call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),&
-                      cons_p(8),cons_p(6),cons_p(7),&
-                      tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5),&
-                      tmp_flux(8),tmp_flux(6),tmp_flux(7), &
-                      vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, &
-                      avg_det,avg_alp,avg_beta)
-                 
-                 densflux(i,j,k) = 0.5d0 * tmp_flux(1)
-                 sxflux(i,j,k)   = 0.5d0 * tmp_flux(2)
-                 syflux(i,j,k)   = 0.5d0 * tmp_flux(3)
-                 szflux(i,j,k)   = 0.5d0 * tmp_flux(4)
-                 tauflux(i,j,k)  = 0.5d0 * tmp_flux(5)
-                 Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6)
-                 Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
-                 Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
-                 
-                 if(clean_divergence.ne.0) then
-                    psidcf = cons_p(8)
-                    psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
-                 endif
-
-                 call num_x_fluxM(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
-                      cons_m(8),cons_m(6),cons_m(7),&
-                      tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5), &
-                      tmp_flux(8),tmp_flux(6),tmp_flux(7), &
-                      vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, &
-                      avg_det,avg_alp,avg_beta)
-                 
-                 densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1)
-                 sxflux(i,j,k)   = sxflux(i,j,k)   + 0.5d0 * tmp_flux(2)
-                 syflux(i,j,k)   = syflux(i,j,k)   + 0.5d0 * tmp_flux(3)
-                 szflux(i,j,k)   = szflux(i,j,k)   + 0.5d0 * tmp_flux(4)
-                 tauflux(i,j,k)  = tauflux(i,j,k)  + 0.5d0 * tmp_flux(5)
-                 Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6)
-                 Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
-                 Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
-                 
-                 if(clean_divergence.ne.0) then
-                    psidcf = cons_m(8)
-                    psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
-                 endif
-
-              else
-                 
-                 call CCTK_WARN(0, "Flux direction not x,y,z")
-                 
-              end if
-              
-           end do
-        end do
-     end do
-     
-     if (CCTK_EQUALS(riemann_solver,"Roe")) then   
-
-        call CCTK_WARN(0, "Roe and Marquina not implemented in MHD yet!!!")
-        
-!!$        
-!!$        call GRHydro_RoeSolveGeneral(CCTK_PASS_FTOF)
-!!$        
-!!$        if (evolve_tracer .ne. 0) then
-!!$           
-!!$           call GRHydro_HLLE_TracerGeneral(CCTK_PASS_FTOF)
-!!$           
-!!$        end if
-!!$        
-     else if (CCTK_EQUALS(riemann_solver,"Marquina")) then   
-
-        call CCTK_WARN(0, "Roe and Marquina not implemented in MHD yet!!!")
-
-!!$        
-!!$        call GRHydro_MarquinaGeneral(CCTK_PASS_FTOF)
-!!$        
-!!$    Tracers are built directly in to the Marquina solver
-        
-     end if
-     
-  end if
-  
-end subroutine RiemannSolveGeneralM