make EOS change work with MP

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

1 files changed, 53 insertions, 3 deletions

diff --git a/src/GRHydro_EoSChangeGamma.F90 b/src/GRHydro_EoSChangeGamma.F90
index 2820c1b..fb075e2 100644
--- a/src/GRHydro_EoSChangeGamma.F90
+++ b/src/GRHydro_EoSChangeGamma.F90
@@ -148,6 +148,12 @@ subroutine GRHydro_EoSChangeK(CCTK_ARGUMENTS)
 
   implicit none
 
+  ! save memory when MP is not used
+  ! TARGET has to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1
+  TARGET gaa, gab, gac, gbb, gbc, gcc
+  TARGET gxx, gxy, gxz, gyy, gyz, gzz
+  TARGET lvel, vel
+
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
 
@@ -158,6 +164,11 @@ subroutine GRHydro_EoSChangeK(CCTK_ARGUMENTS)
   
   CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: Q
 
+  ! save memory when MP is not used
+  CCTK_INT :: GRHydro_UseGeneralCoordinates
+  CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33
+  CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup
+
 !!$  Set up the fluid constants
 ! begin EOS Omni vars
   integer :: n,keytemp,anyerr,keyerr(1)
@@ -174,12 +185,16 @@ subroutine GRHydro_EoSChangeK(CCTK_ARGUMENTS)
   local_Gamma = 1.0d0 + xpress/xeps
   local_K = xpress
 
+  call CCTK_INFO("Adjusting EOS via K change!")
+
   if (abs(local_Gamma - 2.d0) < 1.d-10) then
 
     rho = -0.5d0/local_k+sqrt(0.25d0/local_k**2+(rho+initial_k*rho**2)/local_k)
 
   else if (abs(local_Gamma - 3.d0) < 1.d-10) then
 
+     call CCTK_WARN(0, "This part of the code is not tested!")
+
 !!$ This code is probably just wrong. We have never used it anyway.
     Q = -9.d0 * local_k**2 * rho * (2.d0 + 3.d0 * initial_k * rho**2) + &
          sqrt(local_k**3 * (32.d0 + 81.d0 * local_k * rho**2 * &
@@ -193,11 +208,46 @@ subroutine GRHydro_EoSChangeK(CCTK_ARGUMENTS)
   else
     call CCTK_WARN(0, "EoSChangeK only knows how to do Gamma=2 or 3!")
   end if
-  
+ 
   press = local_k * rho**local_gamma
-  eps = local_k / (local_gamma - 1.d0) * rho**(local_gamma-1)
+  eps = local_k / (local_gamma - 1.d0) * rho**(local_gamma-1.0d0)
+
+  ! save memory when MP is not used
+  if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
+    g11 => gaa
+    g12 => gab
+    g13 => gac
+    g22 => gbb
+    g23 => gbc
+    g33 => gcc
+    vup => lvel
+  else
+    g11 => gxx
+    g12 => gxy
+    g13 => gxz
+    g22 => gyy
+    g23 => gyz
+    g33 => gzz
+    vup => vel
+  end if
+
+  do k = 1, cctk_lsh(3)
+    do j = 1, cctk_lsh(2)
+      do i = 1, cctk_lsh(1)
+        
+        det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),\
+             g22(i,j,k),g23(i,j,k),g33(i,j,k))
+        call prim2conpolytype(GRHydro_polytrope_handle,g11(i,j,k),g12(i,j,k),&
+             g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(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 do
+    end do
+  end do
+
 
-  call Primitive2ConservativePolyCells(CCTK_ARGUMENTS)
 
 end subroutine GRHydro_EoSChangeK