file/parameter string replacement from whisky to GRHydro

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

1 files changed, 372 insertions, 0 deletions

diff --git a/src/GRHydro_RiemannSolve.F90 b/src/GRHydro_RiemannSolve.F90
new file mode 100644
index 0000000..76739a3
--- /dev/null
+++ b/src/GRHydro_RiemannSolve.F90
@@ -0,0 +1,372 @@
+ /*@@
+   @file      GRHydro_RiemannSolve.F90
+   @date      Sat Jan 26 02:20:25 2002
+   @author    
+   @desc 
+   A wrapper routine to call the correct Riemann solver
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+
+ /*@@
+   @routine    RiemannSolve
+   @date       Sat Jan 26 02:20:48 2002
+   @author     Pedro Montero, Ian Hawke
+   @desc 
+   A wrapper routine to switch between the different Riemann solvers.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine RiemannSolve(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i,j,k
+  
+  if (CCTK_EQUALS(riemann_solver,"HLLE")) then
+    
+    call GRHydro_HLLE(CCTK_PASS_FTOF)
+
+    if (evolve_tracer .ne. 0) then
+    
+      call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF)
+
+    end if
+    
+  else if (CCTK_EQUALS(riemann_solver,"Roe")) then   
+    
+    call GRHydro_RoeSolve(CCTK_PASS_FTOF)
+
+    if (evolve_tracer .ne. 0) then
+    
+      call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF)
+
+    end if
+    
+  else if (CCTK_EQUALS(riemann_solver,"Marquina")) then   
+    
+    call GRHydro_Marquina(CCTK_PASS_FTOF)
+
+!!$    Tracers are built directly in to the Marquina solver
+
+  end if
+
+end subroutine RiemannSolve
+
+ /*@@
+   @routine    RiemannSolvePolytype
+   @date       Tue Mar 19 11:40:20 2002
+   @author     Ian Hawke
+   @desc 
+   The same as above, just specializing to polytropic type EOS.
+   Currently there is no point to this routine right now.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+
+subroutine RiemannSolvePolytype(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i,j,k
+  
+  if (CCTK_EQUALS(riemann_solver,"HLLE")) then
+    
+    call GRHydro_HLLE(CCTK_PASS_FTOF)
+
+    if (evolve_tracer .ne. 0) then
+    
+      call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF)
+
+    end if
+    
+  else if (CCTK_EQUALS(riemann_solver,"Roe")) then   
+    
+    call GRHydro_RoeSolve(CCTK_PASS_FTOF)
+
+    if (evolve_tracer .ne. 0) then
+    
+      call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF)
+
+    end if
+    
+  else if (CCTK_EQUALS(riemann_solver,"Marquina")) then   
+    
+    call GRHydro_Marquina(CCTK_PASS_FTOF)
+
+!!$    Tracers are built directly in to the Marquina solver
+
+  end if
+
+end subroutine RiemannSolvePolytype
+
+
+
+ /*@@
+   @routine    RiemannSolveGeneral
+   @date       Tue Mar 19 11:40:20 2002
+   @author     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 RiemannSolveGeneral(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(5) :: tmp_flux, cons_p, cons_m
+  CCTK_REAL, dimension(6) :: prim_p, prim_m
+  CCTK_REAL :: avg_det, avg_alp, avg_beta
+  CCTK_REAL :: psi4h, gxxh, gyyh, gzzh, gxyh, gxzh, gyzh
+
+  densflux = 0.d0
+  sxflux   = 0.d0
+  syflux   = 0.d0
+  szflux   = 0.d0
+  tauflux  = 0.d0
+
+!!$  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_HLLEGeneral(CCTK_PASS_FTOF)
+    
+    if (evolve_tracer .ne. 0) then
+    
+      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_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) 
+          
+          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)
+          
+          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) 
+        
+!!$        Set metric terms at interface
+          
+          if (shift_state .ne. 0) then
+            if (flux_direction == 1) then
+              avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+                   betax(i,j,k))
+            else if (flux_direction == 2) then
+              avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+                   betay(i,j,k))
+            else if (flux_direction == 3) then
+              avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+                   betaz(i,j,k))
+            else
+              call CCTK_WARN(0, "Flux direction not x,y,z")
+            end if
+          else
+            avg_beta = 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))
+          
+          if (conformal_state > 0) then
+            
+            psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4
+            gxxh = gxxh * psi4h
+            gxyh = gxyh * psi4h
+            gxzh = gxzh * psi4h
+            gyyh = gyyh * psi4h
+            gyzh = gyzh * psi4h
+            gzzh = gzzh * psi4h
+            
+          end if
+          
+          call SpatialDeterminant(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,avg_det)
+          
+          if (flux_direction == 1) then
+            
+            call num_x_flux(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),&
+                 tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5), &
+                 prim_p(2),prim_p(6), &
+                 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)
+            
+            call num_x_flux(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
+                 tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5), &
+                 prim_m(2),prim_m(6), &
+                 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)
+            
+          else if (flux_direction == 2) then
+            
+            call num_x_flux(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),&
+                 tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5), &
+                 prim_p(3),prim_p(6), &
+                 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)
+            
+            call num_x_flux(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
+                 tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5), &
+                 prim_m(3),prim_m(6), &
+                 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)
+            
+          else if (flux_direction == 3) then
+            
+            call num_x_flux(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),&
+                 tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5), &
+                 prim_p(4),prim_p(6), &
+                 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)
+            
+            call num_x_flux(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
+                 tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5), &
+                 prim_m(4),prim_m(6), &
+                 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)
+            
+          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 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 GRHydro_MarquinaGeneral(CCTK_PASS_FTOF)
+      
+!!$    Tracers are built directly in to the Marquina solver
+
+    end if
+
+ end if
+
+end subroutine RiemannSolveGeneral
+
+