This is a _temporary_ repository to be able to start to work on the

code right now. I have put in the public version of Whisky to start from.

Everybody with commit rights should get commit messages (and the other
way around). It should not be a problem to add people to that list, just
ask. I don't want to get into political problems because someone feels
excluded, but I also don't want to give everyone access per se.

Frank



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

1 files changed, 174 insertions, 0 deletions

diff --git a/src/Utils.F90 b/src/Utils.F90
new file mode 100644
index 0000000..935a607
--- /dev/null
+++ b/src/Utils.F90
@@ -0,0 +1,174 @@
+ /*@@
+   @file      Utils.F
+   @date      Sat Jan 26 02:28:46 2002
+   @author    
+   @desc 
+   Utility functions for other thorns. Calculation of the determinant
+   of the spatial metric and the upper metric.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+ /*@@
+   @routine    Whisky_RefinementLevel
+   @date       July 2005
+   @author     
+   @desc 
+   Calculates the current refinement level from flesh data
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+@@*/
+
+subroutine Whisky_RefinementLevel(CCTK_ARGUMENTS)
+
+  implicit none
+  DECLARE_CCTK_ARGUMENTS
+
+  whisky_reflevel = aint(log10(dble(cctk_levfac(1)))/log10(2.0d0))
+
+end subroutine Whisky_RefinementLevel
+
+
+ /*@@
+   @routine    SpatialDeterminant
+   @date       Sat Jan 26 02:30:23 2002
+   @author     
+   @desc 
+   Calculates the determinant of the spatial metric.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+@@*/
+
+subroutine SpatialDeterminant(gxx,gxy,gxz,gyy,gyz,gzz,det)
+
+  implicit none
+  
+  CCTK_REAL :: det
+  CCTK_REAL :: gxx,gxy,gxz,gyy,gyz,gzz
+
+  det = -gxz**2*gyy + 2*gxy*gxz*gyz - gxx*gyz**2 - gxy**2*gzz + gxx*gyy*gzz
+
+  
+!!$  Why is this weird order necessary? Search me. It just seemed 
+!!$  to make a really odd NaN go away.
+
+!!$  det =  -gxz**2*gyy + 2*gxy*gxz*gyz 
+!!$  det = det - gxx*gyz**2 - gxy**2*gzz
+!!$  det = det + gxx*gyy*gzz
+  
+  return
+  
+end subroutine SpatialDeterminant
+
+
+
+/*@@
+@routine    UpperMetric
+@date       Sat Jan 26 02:32:26 2002
+@author     
+@desc 
+Calculates the upper metric. The determinant is given, not
+calculated.
+@enddesc 
+@calls     
+@calledby   
+@history 
+
+@endhistory 
+@@*/
+
+subroutine UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, &
+     det, gxx, gxy, gxz, gyy, gyz, gzz)
+  
+  implicit none
+  
+  CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, &
+       gxx, gxy, gxz, gyy, gyz, gzz
+  
+  uxx=(-gyz**2 + gyy*gzz)/det
+  uxy=(gxz*gyz - gxy*gzz)/det
+  uyy=(-gxz**2 + gxx*gzz)/det
+  uxz=(-gxz*gyy + gxy*gyz)/det
+  uyz=(gxy*gxz - gxx*gyz)/det
+  uzz=(-gxy**2 + gxx*gyy)/det
+  
+  return
+  
+end subroutine UpperMetric
+
+ /*@@
+   @routine    SetMetricTemps
+   @date       Wed Mar 23 10:01:08 2005
+   @author     Ian Hawke
+   @desc 
+   If using the new EOS store the metric determinant and inverse metric
+   to avoid repeated computation
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine SetMetricTemps(CCTK_ARGUMENTS)
+
+  implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT :: i, j, k, nx, ny, nz
+  CCTK_REAL :: psi4pt
+
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+
+  do k = 1, nz 
+    do j = 1, ny 
+      do i = 1, nx
+        
+        if (conformal_state .eq. 0) then
+          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),&
+               Whisky_det(i,j,k))
+          call UpperMetric(&
+               Whisky_uxx(i,j,k),Whisky_uxy(i,j,k),Whisky_uxz(i,j,k),&
+               Whisky_uyy(i,j,k),Whisky_uyz(i,j,k),Whisky_uzz(i,j,k),&
+               Whisky_det(i,j,k),&
+               gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
+               gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))        
+        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),&
+               Whisky_det(i,j,k))
+          call UpperMetric(&
+               Whisky_uxx(i,j,k),Whisky_uxy(i,j,k),Whisky_uxz(i,j,k),&
+               Whisky_uyy(i,j,k),Whisky_uyz(i,j,k),Whisky_uzz(i,j,k),&
+               Whisky_det(i,j,k),&
+               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))
+        end if
+
+      end do
+    end do
+  end do
+  
+end subroutine SetMetricTemps
+