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

diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90
new file mode 100644
index 0000000..3753797
--- /dev/null
+++ b/src/GRHydro_UpdateMask.F90
@@ -0,0 +1,247 @@
+ /*@@
+   @file      GRHydro_UpdateMask.F90
+   @date      Wed Mar 13 14:18:38 2002
+   @author    
+   @desc 
+   Alter the update terms if inside the atmosphere or excision region
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.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 GRHydroUpdateAtmosphereMask(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT :: i,j,k
+  CCTK_REAL :: frac
+
+  CCTK_INT :: type_bits, atmosphere
+
+  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
+  call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
+
+  frac = CCTK_DELTA_TIME
+
+  !$OMP PARALLEL DO PRIVATE(j,i)
+  do k = 1, cctk_lsh(3)
+    do j = 1, cctk_lsh(2)
+      do i = 1, cctk_lsh(1)
+        if ( GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0) .or. &
+             (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \
+                                          type_bits, atmosphere)) .or. &
+             (tau(i,j,k) + frac * taurhs(i,j,k) .le. 0.d0) .or. &
+             (dens(i,j,k) + frac * densrhs(i,j,k) .le. 0.d0) ) then
+          densrhs(i,j,k) = 0.0d0
+          srhs(i,j,k,:)   = 0.0d0
+          taurhs(i,j,k)  = 0.0d0
+          atmosphere_mask(i,j,k) = 1
+
+          SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
+        end if
+      end do
+    end do
+  end do
+  !$OMP END PARALLEL DO
+
+end subroutine GRHydroUpdateAtmosphereMask
+
+ /*@@
+   @routine    GRHydro_SetupMask
+   @date       Thu Jun 20 13:27:28 2002
+   @author     Ian Hawke
+   @desc 
+   Initialize the mask to be zero.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_SetupMask(CCTK_ARGUMENTS)
+
+  implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT :: type_bits, not_atmosphere, i, j, k
+
+  ! Initialize all rhs variables and the mask.
+  ! The former vars need to be initialized since there is
+  ! no rhs computation in CCTK_INITIAL or POSTINITIAL.
+
+  densrhs = 0.0d0
+  taurhs  = 0.0d0 
+  srhs    = 0.0d0
+
+  if (evolve_tracer .ne. 0) then
+     cons_tracerrhs = 0.0d0
+  endif
+
+  atmosphere_mask = 0
+  
+  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
+  call SpaceMask_GetStateBits(not_atmosphere, &
+       &"Hydro_Atmosphere", "not_in_atmosphere")
+
+  do k = 1, cctk_lsh(3)
+    do j = 1, cctk_lsh(2)
+      do i = 1, cctk_lsh(1)
+  
+        SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, not_atmosphere)
+
+      end do
+    end do
+  end do
+
+  call CCTK_INFO("Setting up the atmosphere mask: all points are not_atmosphere")
+  
+end subroutine GRHydro_SetupMask
+
+ /*@@
+   @routine    GRHydro_InitAtmosMask
+   @date       Thu Jun 20 13:27:28 2002
+   @author     Ian Hawke
+   @desc 
+   Initialize the mask based on rho_min. This is used only if wk_atmosphere=yes.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_InitAtmosMask(CCTK_ARGUMENTS)
+
+  implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+
+  CCTK_INT :: type_bits, atmosphere, i, j, k
+
+  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
+  call SpaceMask_GetStateBits(atmosphere, &
+       &"Hydro_Atmosphere", "in_atmosphere")
+
+  do k = 1, cctk_lsh(3)
+    do j = 1, cctk_lsh(2)
+      do i = 1, cctk_lsh(1)
+
+        if (rho(i,j,k) .le. GRHydro_rho_min) then
+  
+          SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
+          atmosphere_mask(i,j,k) = 1
+
+        end if
+        
+      end do
+    end do
+  end do
+
+  call CCTK_INFO("Setting up the atmosphere mask: points with rho<rho_min are set to be atmosphere")
+  
+end subroutine GRHydro_InitAtmosMask
+
+ /*@@
+   @routine    GRHydro_AtmosphereReset
+   @date       Thu Jun 20 13:30:51 2002
+   @author     Ian Hawke
+   @desc 
+   After MoL has evolved, if a point is supposed to be reset then do so.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS)
+
+  implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT :: i, j, k
+  CCTK_REAL :: det, psi4pt
+
+  CCTK_INT :: type_bits, atmosphere, not_atmosphere
+  
+  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
+  call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",&
+                              "in_atmosphere")
+  call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",&
+                              "not_in_atmosphere")
+
+  do k = 1, cctk_lsh(3)
+    do j = 1, cctk_lsh(2)
+      do i = 1, cctk_lsh(1)
+        
+        if ( (atmosphere_mask(i, j, k) .eq. 1) &
+             &.or. (SpaceMask_CheckStateBitsF90(space_mask,i, j, k, type_bits,\
+                    atmosphere)) &
+             &) then
+
+!!$          write(*,*) 'Resetting at ',i,j,k, atmosphere_mask(i, j, k), &
+!!$             & (SpaceMask_CheckStateBitsF90(space_mask,i, j, k, type_bits,\
+!!$                    atmosphere)) 
+
+          rho(i,j,k) = GRHydro_rho_min
+          velx(i,j,k) = 0.0d0
+          vely(i,j,k) = 0.0d0
+          velz(i,j,k) = 0.0d0
+          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), det)
+            call prim2conpolytype(GRHydro_polytrope_handle, &
+                 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, &
+                 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), 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))
+          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)
+            call prim2conpolytype(GRHydro_polytrope_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), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
+                 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
+          if (wk_atmosphere .eq. 0) then
+            atmosphere_mask(i, j, k) = 0
+            SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\
+                                      not_atmosphere)
+          end if
+
+        end if
+
+      end do
+    end do
+  end do
+
+!!$  call GRHydro_Boundaries(CCTK_PASS_FTOF)
+
+end subroutine GRHydro_AtmosphereReset