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

diff --git a/src/GRHydro_Particle.F90 b/src/GRHydro_Particle.F90
new file mode 100644
index 0000000..464bf78
--- /dev/null
+++ b/src/GRHydro_Particle.F90
@@ -0,0 +1,248 @@
+ /*@@
+   @file      GRHydro_Particle.F90
+   @date      Wed Mar 31 11:10:52 2004
+   @author    Ian Hawke
+   @desc 
+   Track coordinates of particles.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+
+ /*@@
+   @routine    GRHydroParticleRHS
+   @date       Wed Mar 31 11:11:26 2004
+   @author     Ian Hawke
+   @desc 
+   Track particles.
+   Most of this code is taken from Peter Dieners EHFinder - generator tracking.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydroParticleRHS(CCTK_ARGUMENTS)
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i, l
+  CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
+
+  character(len=200) :: particle_interp
+  character(len=128) :: warn_message
+  CCTK_INT :: particle_interp_len
+  character(len=7) :: particle_order
+  character(len=15) :: vname
+
+  CCTK_INT, dimension(1) :: lsh
+  CCTK_POINTER, dimension(3) :: interp_coords
+  CCTK_POINTER, dimension(7) :: out_arrays
+  CCTK_INT, dimension(7) :: in_arrays
+  CCTK_INT, dimension(7), parameter :: op_indices = (/ 0, 1, 2, 3, 4, &
+                                                        5, 6 /), &
+                                        op_codes = (/ 0, 0, 0, 0, 0, &
+                                                      0, 0 /)
+  CCTK_INT, dimension(7) :: out_types
+
+  out_types = CCTK_VARIABLE_REAL
+
+! Convert the particle_interpolator string parameter to a Fortran string.
+  call CCTK_FortranString ( particle_interp_len, particle_interpolator, &
+       particle_interp )
+
+! Get the corresponding interpolator handle.
+  call CCTK_InterpHandle ( interp_handle, particle_interp )
+
+  if ( interp_handle .lt. 0 ) then
+    warn_message = 'Cannot get handle for interpolation.'
+    warn_message = trim(warn_message)//' Forgot to activate an implementation'
+    warn_message = trim(warn_message)//' providing interpolation operators?'
+    call CCTK_WARN( 0, trim(warn_message) )
+  end if
+
+! Convert the interpolation order parameter to a Fortran string to be placed
+! in the interpolator table. Note that the order is assumed to contain only
+! 1 digit.
+  write(particle_order,'(a6,i1)') 'order=', particle_interpolation_order
+
+! Create the table directly from the string.
+  call Util_TableCreateFromString ( table_handle, particle_order )
+  if ( table_handle .lt. 0 ) then
+    call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' )
+  end if
+
+! Get the 3D coordinate system handle.
+  call CCTK_CoordSystemHandle ( coord_system_handle, 'cart3d' )
+  if ( coord_system_handle .lt. 0) then
+    warn_message = 'Cannot get handle for cart3d coordinate system.'
+    warn_message = trim(warn_message)//' Forgot to activate an implementation'
+    warn_message = trim(warn_message)//' providing coordinates?'
+    call CCTK_WARN( 0, trim(warn_message) )
+  endif
+
+! Find out how many interpolation points are located on this processor.
+  call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, 'GRHydro::particles' )
+  if ( status .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get local size for particles' )
+  end if
+
+! Set the pointers to the output arrays.
+  out_arrays(1) = CCTK_PointerTo(particle_vx)
+  out_arrays(2) = CCTK_PointerTo(particle_vy)
+  out_arrays(3) = CCTK_PointerTo(particle_vz)
+  out_arrays(4) = CCTK_PointerTo(particle_alp)
+  out_arrays(5) = CCTK_PointerTo(particle_betax)
+  out_arrays(6) = CCTK_PointerTo(particle_betay)
+  out_arrays(7) = CCTK_PointerTo(particle_betaz)
+
+!   Set the pointers to the points to be interpolated to.
+  interp_coords(1) = CCTK_PointerTo(particle_x)
+  interp_coords(2) = CCTK_PointerTo(particle_y)
+  interp_coords(3) = CCTK_PointerTo(particle_z)
+
+!   Set the indices to the input grid functions.
+  call CCTK_VarIndex ( in_arrays(1), 'GRHydro::velx' )
+  call CCTK_VarIndex ( in_arrays(2), 'GRHydro::vely' )
+  call CCTK_VarIndex ( in_arrays(3), 'GRHydro::velz' )
+  call CCTK_VarIndex ( in_arrays(4), 'admbase::alp' )
+  call CCTK_VarIndex ( in_arrays(5), 'admbase::betax' )
+  call CCTK_VarIndex ( in_arrays(6), 'admbase::betay' )
+  call CCTK_VarIndex ( in_arrays(7), 'admbase::betaz' )
+
+!     Set the operand indices table entry
+  call Util_TableSetIntArray ( status, table_handle, 7, &
+       op_indices(1:7), 'operand_indices' )
+  if ( status .lt. 0 ) then
+    warn_message = 'Cannot set operand indices array in parameter table' 
+    call CCTK_WARN ( 0, trim(warn_message) )
+  endif
+
+!     Set the corresponding table entry for the operation codes.
+  call Util_TableSetIntArray ( status, table_handle, 7, &
+       op_codes(1:7), 'operation_codes' )
+  if ( status .lt. 0 ) then
+    warn_message = 'Cannot set operation codes array in parameter table'
+    call CCTK_WARN ( 0, trim(warn_message) )
+  endif
+
+!     Call the interpolator.
+  call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
+       table_handle, coord_system_handle, &
+       lsh(1), CCTK_VARIABLE_REAL, &
+       interp_coords, 7, in_arrays(1:7), &
+       7, out_types(1:7), out_arrays(1:7) )
+
+  if ( status .lt. 0 ) then
+    call CCTK_INFO ( 'Interpolation failed.' )
+  end if
+
+!     For each point on this processor calculate the right hand side of the
+!     characteristic evolution equation.
+  do i = 1, lsh(1)
+
+    particle_x_rhs(i) = particle_alp(i) * particle_vx(i) - &
+         particle_betax(i)
+    particle_y_rhs(i) = particle_alp(i) * particle_vy(i) - &
+         particle_betay(i)
+    particle_z_rhs(i) = particle_alp(i) * particle_vz(i) - &
+         particle_betaz(i)
+
+  end do
+  
+  return
+
+end subroutine GRHydroParticleRHS
+
+/*@@
+   @routine    GRHydroParticleInitial
+   @date       Wed Mar 31 11:53:25 2004
+   @author     Ian Hawke
+   @desc 
+   Sets up the initial coordinates of the particles. Can be overwritten
+   by something scheduled in HydroBase_Initial.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydroParticleInitial(CCTK_ARGUMENTS)
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT, dimension(1) :: lsh, lbnd
+  CCTK_INT :: i, status
+  CCTK_REAL :: radius, angle, xmin, xmax, ymin, ymax, zmin, zmax, twopi
+
+  twopi = 4.d0 * atan(1.d0)
+
+  if (number_of_particles .gt. 0) then
+
+! Find out how many interpolation points are located on this processor.
+    call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, 'GRHydro::particles' )
+    if ( status .lt. 0 ) then
+      call CCTK_WARN ( 0, 'cannot get local size for particles' )
+    end if
+    call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, 'GRHydro::particles' )
+    if ( status .lt. 0 ) then
+      call CCTK_WARN ( 0, 'cannot get lower bounds for particles' )
+    end if
+
+    call CCTK_CoordRange(status,cctkGH,xmin,xmax,-1,"x","cart3d")
+    call CCTK_CoordRange(status,cctkGH,ymin,ymax,-1,"y","cart3d")
+    call CCTK_CoordRange(status,cctkGH,zmin,zmax,-1,"z","cart3d")
+
+
+    radius = 0.1d0 * min(xmax, ymax, zmax)
+
+    do i = 1, lsh(1)
+
+      angle = dble(lbnd(1)+i-1) / dble(number_of_particles) * twopi 
+
+      particle_x(i) = radius * sin(angle)
+      particle_y(i) = radius * cos(angle)
+      particle_z(i) = 0.d0
+
+      particle_x_p(i) = particle_x(i)
+      particle_y_p(i) = particle_y(i)
+      particle_z_p(i) = particle_z(i)
+
+      particle_x_p_p(i) = particle_x(i)
+      particle_y_p_p(i) = particle_y(i)
+      particle_z_p_p(i) = particle_z(i)
+
+      particle_x_rhs(i) = 0.d0
+      particle_y_rhs(i) = 0.d0
+      particle_z_rhs(i) = 0.d0
+      
+      particle_vx(i) = 0.d0
+      particle_vy(i) = 0.d0
+      particle_vz(i) = 0.d0
+      particle_alp(i) = 0.d0
+      particle_betax(i) = 0.d0
+      particle_betay(i) = 0.d0
+      particle_betaz(i) = 0.d0
+
+    end do
+    
+  end if
+  
+end subroutine GRHydroParticleInitial