add Multipatch support to GRHydro

* not all features of GRHydro are supported yet, in particular only the
  HLLE solver supports Mulitpatch yet.

Original commit by Christian Reisswig and Christian Ott



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

1 files changed, 57 insertions, 57 deletions

diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90
index d3c9732..3c51302 100644
--- a/src/GRHydro_Con2Prim.F90
+++ b/src/GRHydro_Con2Prim.F90
@@ -101,11 +101,11 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
          
         epsnegative = .false.
 
-        det = SPATIAL_DETERMINANT(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 = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k), \
+              gbb(i,j,k),gbc(i,j,k),gcc(i,j,k))
         call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,&
-             gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
-             gyz(i,j,k),gzz(i,j,k))        
+             gaa(i,j,k),gab(i,j,k),gac(i,j,k),gbb(i,j,k),&
+             gbc(i,j,k),gcc(i,j,k))        
 
 !!$        if (det < 0.e0) then  
 !!$          call CCTK_WARN(0, "nan produced (det negative)")
@@ -135,7 +135,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
            dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
            rho(i,j,k) = GRHydro_rho_min
            scon(i,j,k,:) = 0.d0
-           vel(i,j,k,:) = 0.d0
+           lvel(i,j,k,:) = 0.d0
            w_lorentz(i,j,k) = 1.d0
 
            if(evolve_temper.ne.0) then
@@ -148,7 +148,6 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
                    rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
                    press(i,j,k),keyerr,anyerr)
            else
-              ! w_lorentz=1, so the expression for tau reduces to:
               keytemp = 0
               call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                    rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
@@ -157,6 +156,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
                    rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
            endif
 
+           ! w_lorentz=1, so the expression for tau reduces to:
            tau(i,j,k)  = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k)
            
            cycle
@@ -165,15 +165,15 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
 
          if(evolve_temper.eq.0) then
             call Con2Prim_pt(GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2), &
-                 scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vel(i,j,k,1),vel(i,j,k,2), &
-                 vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
+                 scon(i,j,k,3),tau(i,j,k),rho(i,j,k),lvel(i,j,k,1),lvel(i,j,k,2), &
+                 lvel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
                  uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
                  z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & 
                  GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
          else
             call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),&
-                 scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),vel(i,j,k,1),&
-                 vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
+                 scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),lvel(i,j,k,1),&
+                 lvel(i,j,k,2),lvel(i,j,k,3),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
                  press(i,j,k),w_lorentz(i,j,k), &
                  uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
                  z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & 
@@ -186,7 +186,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
                local_perc_ptol = GRHydro_perc_ptol*100.0d0
                call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),&
                     scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),&
-                    vel(i,j,k,1),vel(i,j,k,2), vel(i,j,k,3),eps(i,j,k),&
+                    lvel(i,j,k,1),lvel(i,j,k,2), lvel(i,j,k,3),eps(i,j,k),&
                     temperature(i,j,k),y_e(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
                     uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
                     z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & 
@@ -237,7 +237,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
           dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
           rho(i,j,k) = GRHydro_rho_min
           scon(i,j,k,:) = 0.d0
-          vel(i,j,k,:) = 0.d0
+          lvel(i,j,k,:) = 0.d0
           w_lorentz(i,j,k) = 1.d0
 
            call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
@@ -259,8 +259,8 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
           !$OMP END CRITICAL
           call Con2Prim_ptPolytype(GRHydro_polytrope_handle, &
                dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
-               tau(i,j,k), rho(i,j,k), vel(i,j,k,1), vel(i,j,k,2), &
-               vel(i,j,k,3), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k), &
+               tau(i,j,k), rho(i,j,k), lvel(i,j,k,1), lvel(i,j,k,2), &
+               lvel(i,j,k,3), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k), &
                uxx, uxy, uxz, uyy, uyz, uzz, det, x(i,j,k), y(i,j,k), &
                z(i,j,k), r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
 #endif          
@@ -1048,18 +1048,18 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS)
         if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
             GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
          
-        gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
-        gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
-        gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
-        gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
-        gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
-        gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
-        gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
-        gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
-        gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
-        gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
-        gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
-        gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+        gxxl = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset))
+        gxyl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset))
+        gxzl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset))
+        gyyl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset))
+        gyzl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset))
+        gzzl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset))
+        gxxr = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset))
+        gxyr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset))
+        gxzr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset))
+        gyyr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset))
+        gyzr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset))
+        gzzr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset))
 
         epsnegative = .false.
 
@@ -1251,11 +1251,11 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS)
         if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
             GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
 
-        det = SPATIAL_DETERMINANT(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 = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k),\
+             gbb(i,j,k),gbc(i,j,k),gcc(i,j,k))
         call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,&
-             gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
-             gyz(i,j,k),gzz(i,j,k))        
+             gaa(i,j,k),gab(i,j,k),gac(i,j,k),gbb(i,j,k),&
+             gbc(i,j,k),gcc(i,j,k))        
 
         if (evolve_tracer .ne. 0) then
            do itracer=1,number_of_tracers
@@ -1277,8 +1277,8 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS)
 
         call Con2Prim_ptPolytype(GRHydro_eos_handle, dens(i,j,k),&
              scon(i,j,k,1),scon(i,j,k,2), &
-             scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vel(i,j,k,1),vel(i,j,k,2), &
-             vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
+             scon(i,j,k,3),tau(i,j,k),rho(i,j,k),lvel(i,j,k,1),lvel(i,j,k,2), &
+             lvel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
              uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
              z(i,j,k),r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
                 
@@ -1600,18 +1600,18 @@ subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS)
         if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
             GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
         
-        gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
-        gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
-        gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
-        gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
-        gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
-        gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
-        gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
-        gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
-        gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
-        gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
-        gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
-        gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+        gxxl = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset))
+        gxyl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset))
+        gxzl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset))
+        gyyl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset))
+        gyzl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset))
+        gzzl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset))
+        gxxr = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset))
+        gxyr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset))
+        gxzr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset))
+        gyyr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset))
+        gyzr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset))
+        gzzr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset))
 
         avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
         avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
@@ -1686,18 +1686,18 @@ subroutine Con2PrimBoundsTracer(CCTK_ARGUMENTS)
     do j = GRHydro_stencil, ny - GRHydro_stencil + 1
       do i = GRHydro_stencil, nx - GRHydro_stencil + 1
         
-        gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
-        gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
-        gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
-        gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
-        gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
-        gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
-        gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
-        gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
-        gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
-        gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
-        gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
-        gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+        gxxl = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset))
+        gxyl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset))
+        gxzl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset))
+        gyyl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset))
+        gyzl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset))
+        gzzl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset))
+        gxxr = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset))
+        gxyr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset))
+        gxzr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset))
+        gyyr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset))
+        gyzr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset))
+        gzzr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset))
 
         avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
         avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
@@ -1878,7 +1878,7 @@ subroutine check_GRHydro_C2P_failed(CCTK_ARGUMENTS)
           write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k)
           call CCTK_WARN(1,warnline)
           write(warnline,'(a20,3g16.7)') 'velocities: ',&
-               vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3)
+               lvel(i,j,k,1),lvel(i,j,k,2),lvel(i,j,k,3)
           call CCTK_WARN(0,warnline)
           !$OMP END CRITICAL