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

diff --git a/src/GRHydro_Tmunu.F90 b/src/GRHydro_Tmunu.F90
new file mode 100644
index 0000000..6ba34e3
--- /dev/null
+++ b/src/GRHydro_Tmunu.F90
@@ -0,0 +1,175 @@
+ /*@@
+   @file      GRHydro_Tmunu.F90
+   @date      Thu Apr 16 19:38:40 2009
+   @author    Ian Hawke
+   @histpry
+    Apr. 2009: Luca Baiotti copied and adapted for the Tmunu-thorn mechanism the original include file
+   @desc 
+     The calculation of the stress energy tensor.
+     The version used here was worked out by Miguel Alcubierre. I
+     think it was an extension of the routine from GR3D, written
+     by Mark Miller.
+     C version added by Ian Hawke.
+
+     Lower components of the stress-energy tensor obtained from
+     the hydro variables.  The components are given by:
+
+     T      =  rho h  u   u    +  P  g
+      mu nu            mu  nu         mu nu
+
+     where rho is the energy density of the fluid, h the enthalpy
+     and P the pressure.  The enthalpy is given in terms of the
+     basic variables as:
+
+     h  =  1  +  e  +  P/rho
+
+     with e the internal energy (eps here).
+
+     In the expresion for T_{mu,nu} we also have the four-velocity
+     of the fluid given by (v_i is the 3-velocity field):
+
+                                 i
+     u  =  W ( - alpha +  v  beta  )
+      0                    i
+
+     u  =  W v
+      i       i
+                                                i  -1/2
+     with W the Lorentz factor:   W = ( 1 -  v v  )
+                                              i
+
+     and where alpha and beta are the lapse and shift vector.
+
+     Finally, the 4 metric is given by
+
+                   2             i
+     g   =  - alpha  + beta  beta
+      00                   i
+
+     g   =  beta
+      0i        i
+
+
+     g   =  gamma      (the spatial metric)
+      ij        ij
+
+
+   @enddesc 
+ @@*/
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "SpaceMask.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)
+
+ subroutine GRHydro_Tmunu(CCTK_ARGUMENTS)
+
+   implicit none
+
+   DECLARE_CCTK_ARGUMENTS
+
+   CCTK_REAL velxlow, velylow, velzlow
+   CCTK_REAL betaxlow, betaylow, betazlow, beta2
+   CCTK_REAL utlow, uxlow, uylow, uzlow
+   CCTK_REAL rhoenthalpy
+   CCTK_REAL ut,ux,uy,uz,bst,bsx,bsy,bsz,bs2
+   CCTK_INT i,j,k
+
+
+   !$OMP PARALLEL DO PRIVATE(i,j,velxlow, velylow, velzlow,&
+   !$OMP betaxlow, betaylow, betazlow, beta2, utlow, uxlow, uylow, uzlow,&
+   !$OMP rhoenthalpy, ut,ux,uy,uz,bst,bsx,bsy,bsz,bs2)
+   do k = 1, cctk_lsh(3)
+     do j = 1, cctk_lsh(2)
+       do i = 1, cctk_lsh(1)
+
+         if (conformal_state .eq. 0) then
+           velxlow = gxx(i,j,k)*velx(i,j,k) + gxy(i,j,k)*vely(i,j,k) + gxz(i,j,k)*velz(i,j,k)
+           velylow = gxy(i,j,k)*velx(i,j,k) + gyy(i,j,k)*vely(i,j,k) + gyz(i,j,k)*velz(i,j,k)
+           velzlow = gxz(i,j,k)*velx(i,j,k) + gyz(i,j,k)*vely(i,j,k) + gzz(i,j,k)*velz(i,j,k)
+         else
+           velxlow = psi(i,j,k)**4 * (gxx(i,j,k)*velx(i,j,k) + gxy(i,j,k)*vely(i,j,k) + gxz(i,j,k)*velz(i,j,k) )
+           velylow = psi(i,j,k)**4 * (gxy(i,j,k)*velx(i,j,k) + gyy(i,j,k)*vely(i,j,k) + gyz(i,j,k)*velz(i,j,k) )
+           velzlow = psi(i,j,k)**4 * (gxz(i,j,k)*velx(i,j,k) + gyz(i,j,k)*vely(i,j,k) + gzz(i,j,k)*velz(i,j,k) )
+         end if
+
+!!$       Calculate lower components and square of shift vector.
+
+         if (shift_state .ne. 0) then
+
+           if (conformal_state .eq. 0) then
+             betaxlow = gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k)
+             betaylow = gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k)
+             betazlow = gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k)
+           else
+             betaxlow = psi(i,j,k)**4 * (gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k) )
+             betaylow = psi(i,j,k)**4 * (gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k) )
+             betazlow = psi(i,j,k)**4 * (gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k) )
+           end if
+
+           beta2 = betax(i,j,k)*betaxlow + betay(i,j,k)*betaylow + betaz(i,j,k)*betazlow 
+
+         else
+
+           betaxlow = 0.0D0
+           betaylow = 0.0D0
+           betazlow = 0.0D0
+
+           beta2 = 0.0D0
+
+         end if
+
+!!$       Calculate the specific relativistic enthalpy times rho times the
+!!$       square of the lorentz factor.
+
+         rhoenthalpy = w_lorentz(i,j,k)**2*(rho(i,j,k)*(1.0D0 + eps(i,j,k)) + press(i,j,k))
+
+!!$       Calculate lower components of 4-velocity (without the Lorent factor).
+
+         utlow = (-alp(i,j,k) + velx(i,j,k)*betaxlow + vely(i,j,k)*betaylow + velz(i,j,k)*betazlow)
+
+         uxlow = velxlow
+         uylow = velylow
+         uzlow = velzlow
+
+!!$       Calculate Tmunu (the lower components!).
+
+         eTtt(i,j,k) = eTtt(i,j,k) + rhoenthalpy*utlow**2 + press(i,j,k)*(beta2 - alp(i,j,k)**2)
+
+         eTtx(i,j,k) = eTtx(i,j,k) + rhoenthalpy*utlow*uxlow + press(i,j,k)*betaxlow
+         eTty(i,j,k) = eTty(i,j,k) + rhoenthalpy*utlow*uylow + press(i,j,k)*betaylow
+         eTtz(i,j,k) = eTtz(i,j,k) + rhoenthalpy*utlow*uzlow + press(i,j,k)*betazlow
+
+         if (conformal_state .eq. 0) then
+
+           eTxx(i,j,k) = eTxx(i,j,k) + rhoenthalpy*uxlow**2 + press(i,j,k)*gxx(i,j,k)
+           eTyy(i,j,k) = eTyy(i,j,k) + rhoenthalpy*uylow**2 + press(i,j,k)*gyy(i,j,k)
+           eTzz(i,j,k) = eTzz(i,j,k) + rhoenthalpy*uzlow**2 + press(i,j,k)*gzz(i,j,k)
+           
+           eTxy(i,j,k) = eTxy(i,j,k) + rhoenthalpy*uxlow*uylow + press(i,j,k)*gxy(i,j,k)
+           eTxz(i,j,k) = eTxz(i,j,k) + rhoenthalpy*uxlow*uzlow + press(i,j,k)*gxz(i,j,k)
+           eTyz(i,j,k) = eTyz(i,j,k) + rhoenthalpy*uylow*uzlow + press(i,j,k)*gyz(i,j,k)
+
+         else
+
+           eTxx(i,j,k) = eTxx(i,j,k) + rhoenthalpy*uxlow**2 + psi(i,j,k)**4 * press(i,j,k)*gxx(i,j,k)
+           eTyy(i,j,k) = eTyy(i,j,k) + rhoenthalpy*uylow**2 + psi(i,j,k)**4 * press(i,j,k)*gyy(i,j,k)
+           eTzz(i,j,k) = eTzz(i,j,k) + rhoenthalpy*uzlow**2 + psi(i,j,k)**4 * press(i,j,k)*gzz(i,j,k)
+
+           eTxy(i,j,k) = eTxy(i,j,k) + rhoenthalpy*uxlow*uylow + psi(i,j,k)**4 * press(i,j,k)*gxy(i,j,k)
+           eTxz(i,j,k) = eTxz(i,j,k) + rhoenthalpy*uxlow*uzlow + psi(i,j,k)**4 * press(i,j,k)*gxz(i,j,k)
+           eTyz(i,j,k) = eTyz(i,j,k) + rhoenthalpy*uylow*uzlow + psi(i,j,k)**4 * press(i,j,k)*gyz(i,j,k)
+
+         end if
+
+       end do
+     end do
+   end do
+   !$OMP END PARALLEL DO
+
+   return
+ 
+ end subroutine GRHydro_Tmunu