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, 326 insertions, 0 deletions

diff --git a/src/GRHydro_Boundaries.F90 b/src/GRHydro_Boundaries.F90
new file mode 100644
index 0000000..5e8c093
--- /dev/null
+++ b/src/GRHydro_Boundaries.F90
@@ -0,0 +1,326 @@
+ /*@@
+   @file      GRHydro_Boundaries.F90
+   @date      Sat Jan 26 01:01:14 2002
+   @author    
+   @desc 
+   The two routines for dealing with boundary conditions.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+
+#include "util_Table.h"
+
+#define velx(i,j,k) vel(i,j,k,1)
+#define vely(i,j,k) vel(i,j,k,2)
+#define velz(i,j,k) vel(i,j,k,3)
+#define sx(i,j,k) scon(i,j,k,1)
+#define sy(i,j,k) scon(i,j,k,2)
+#define sz(i,j,k) scon(i,j,k,3)
+
+ /*@@
+   @routine    GRHydro_InitSymBound
+   @date       Sat Jan 26 01:03:04 2002
+   @author     Ian Hawke
+   @desc 
+   Sets up the symmetries at the boundaries of the hydrodynamical variables.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   Direct translation of routines from GR3D, GRAstro_Hydro, 
+   written by Mark Miller, or WaveToy routines, or...
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  integer, dimension(3) :: sym
+  integer :: ierr
+  integer :: itracer
+
+  character(len=100) tracername
+  character(len=100) tracerindex
+
+  sym(1) = 1
+  sym(2) = 1
+  sym(3) = 1
+  
+  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::rho")
+  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::press")
+  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::dens")
+  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::tau")
+  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::w_lorentz")
+  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::eps")
+  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::GRHydro_C2P_failed")
+
+!!$ handle multiple tracer variables
+  if(evolve_tracer.eq.1) then
+     call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::GRHydro_tracers")
+     call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::GRHydro_cons_tracers")
+  endif 
+  
+  sym(1) = -1
+  sym(2) = 1
+  sym(3) = 1
+  
+  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[0]")
+  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[0]")
+  
+  sym(1) = 1
+  sym(2) = -1
+  sym(3) = 1
+  
+  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[1]")
+  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[1]")
+  
+  sym(1) = 1
+  sym(2) = 1
+  sym(3) = -1
+  
+  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[2]")
+  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[2]")
+  
+end subroutine GRHydro_InitSymBound
+
+ /*@@
+   @routine    GRHydro_Boundaries
+   @date       Sat Jan 26 01:04:04 2002
+   @author     
+   @desc 
+   Calls the appropriate boundary routines
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_Boundaries(CCTK_ARGUMENTS)
+  
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  integer, dimension(3) :: sw
+  integer :: ierr
+  CCTK_REAL :: pi = 3.141569d0
+  integer :: i,j,k
+
+  CCTK_INT, parameter :: faces=CCTK_ALL_FACES
+  CCTK_INT, parameter :: ione=1
+  
+  sw = GRHydro_stencil
+
+!!$Flat boundaries if required  
+
+  if (CCTK_EQUALS(bound,"flat")) then
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "GRHydro::dens", "Flat")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "GRHydro::tau", "Flat")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "GRHydro::scon", "Flat")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "GRHydro::w_lorentz", "Flat")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "HydroBase::rho", "Flat")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "HydroBase::press", "Flat")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "HydroBase::eps", "Flat")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "HydroBase::vel", "Flat")
+
+    if(evolve_tracer .ne. 0) then 
+      ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+            "GRHydro::GRHydro_tracers", "Flat")
+      ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+            "GRHydro::GRHydro_cons_tracers", "Flat")
+    endif
+
+  endif
+
+  if (CCTK_EQUALS(bound,"none")) then
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "GRHydro::dens", "None")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "GRHydro::tau", "None")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "GRHydro::scon", "None")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "GRHydro::w_lorentz", "None")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "HydroBase::rho", "None")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "HydroBase::press", "None")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "HydroBase::eps", "None")
+    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+         "HydroBase::vel", "None")
+
+    if(evolve_tracer .ne. 0) then 
+      ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+            "GRHydro::GRHydro_tracers", "None")
+      ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
+            "GRHydro::GRHydro_cons_tracers", "None")
+    endif
+
+  end if
+
+  if (CCTK_EQUALS(bound,"scalar")) then
+    call CCTK_WARN(0, "Until somebody uses this I see no reason to support it")
+  end if
+
+  if (ierr < 0) call CCTK_WARN(0, "problems with applying the chosen boundary condition")
+
+end subroutine GRHydro_Boundaries
+
+
+ /*@@
+   @routine    GRHydro_OutflowBoundaries
+   @date       Tue May 17 16:58:06 2005
+   @author     Ian Hawke
+   @desc 
+   Set outflow boundaries over only part of the domain.
+   This is designed to be used with GZPatchSystem.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_OutflowBoundaries(CCTK_ARGUMENTS)
+  
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  integer, dimension(3) :: sw
+  integer :: ierr
+  integer :: i,j,k
+  
+  CCTK_REAL, dimension(3) :: posn
+
+  CCTK_REAL :: det,psi4pt
+
+  sw = GRHydro_stencil
+
+  if (r(1,1,1) < r(1,1,cctk_lsh(3))) then
+
+    if (cctk_bbox(6) .ne. 0) then
+
+      do k = cctk_lsh(3) - sw(3), cctk_lsh(3)
+        do j = 1, cctk_lsh(2)
+          do i = 1, cctk_lsh(1)
+            
+            posn(1) = x(i,j,k)
+            posn(2) = y(i,j,k)
+            posn(3) = z(i,j,k)
+            
+            if (dot_product(outflowboundary_normal, posn) > 0.d0) then
+              
+              rho(i,j,k)        = rho(i,j,k-1)
+              velx(i,j,k)       = velx(i,j,k-1)
+              vely(i,j,k)       = vely(i,j,k-1)
+              velz(i,j,k)       = velz(i,j,k-1)
+              eps(i,j,k)        = eps(i,j,k-1)
+              press(i,j,k)      = press(i,j,k-1)
+              w_lorentz(i,j,k)  = w_lorentz(i,j,k-1)
+              
+              if (conformal_state .eq. 0) then
+                psi4pt = 1d0
+                call SpatialDeterminant(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
+                     gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),det)
+              else
+                psi4pt = psi(i,j,k)**4
+                call SpatialDeterminant(psi4pt*gxx(i,j,k),psi4pt*gxy(i,j,k),&
+                     psi4pt*gxz(i,j,k),psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),&
+                     psi4pt*gzz(i,j,k),det)
+              end if
+              call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),&
+                   psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),&
+                   psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),&
+                   det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+                   tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+                   eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+              
+            end if
+          
+          end do
+        end do
+      end do
+      
+    end if
+  
+  else
+
+    if (cctk_bbox(5) .ne. 0) then
+
+      do k = sw(3), 1, -1
+        do j = 1, cctk_lsh(2)
+          do i = 1, cctk_lsh(1)
+            
+            posn(1) = x(i,j,k)
+            posn(2) = y(i,j,k)
+            posn(3) = z(i,j,k)
+            
+            if (dot_product(outflowboundary_normal, posn) > 0.d0) then
+              
+              rho(i,j,k)        = rho(i,j,k+1)
+              velx(i,j,k)       = velx(i,j,k+1)
+              vely(i,j,k)       = vely(i,j,k+1)
+              velz(i,j,k)       = velz(i,j,k+1)
+              eps(i,j,k)        = eps(i,j,k+1)
+              press(i,j,k)      = press(i,j,k+1)
+              w_lorentz(i,j,k)  = w_lorentz(i,j,k+1)
+              
+              if (conformal_state .eq. 0) then
+                psi4pt = 1d0
+                call SpatialDeterminant(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
+                     gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),det)
+              else
+                psi4pt = psi(i,j,k)**4
+                call SpatialDeterminant(psi4pt*gxx(i,j,k),psi4pt*gxy(i,j,k),&
+                     psi4pt*gxz(i,j,k),psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),&
+                     psi4pt*gzz(i,j,k),det)
+              end if
+              call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),&
+                   psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),&
+                   psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),&
+                   det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+                   tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+                   eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+              
+            end if
+          
+          end do
+        end do
+      end do
+      
+    end if
+
+  end if
+  
+end subroutine GRHydro_OutflowBoundaries
+
+
+