*** empty log message ***

git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/CartGrid3D/trunk@2 c78560ca-4b45-4335-b268-5f3340f3cb52

7 files changed, 480 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..625d155
--- /dev/null
+++ b/README
@@ -0,0 +1,20 @@
+Cactus Code Thorn CartGrid3D
+Authors    : ...
+Managed by : ... <...@...........>
+Version    : ...  
+CVS info   : $Header$
+--------------------------------------------------------------------------
+
+1. Purpose of the thorn
+
+This thorn does ...
+
+2. Dependencies of the thorn
+
+This thorn additionally requires thorns ...
+
+3. Thorn distribution
+
+This thorn is available to ...
+
+4. Additional information
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..100bdb9
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,23 @@
+# Interface definition for thorn CartGrid3D
+# $Header$
+
+implements: grid
+friend: driver
+
+public: 
+ 
+real gridspacings type=SCALAR
+{
+ coarse_dx,
+ coarse_dy,
+ coarse_dz
+} "3D Cartesian grid spacings"
+
+real coordinates type = GF
+{
+ x,
+ y, 
+ z,
+ r
+} "3D Cartesian grid coordinates"
+
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..4832cfa
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,52 @@
+# Parameter definitions for thorn CartGrid3D
+# $Header$
+
+public:
+
+REAL grid_dx "Grid spacing in x-direction"
+{
+  0: ::
+} 0.3
+REAL grid_dy "Grid spacing in y-direction"
+{
+  0: ::
+} 0.3
+REAL grid_dz "Grid spacing in z-direction"
+{
+  0: ::
+} 0.3
+
+REAL grid_xmin "Coordinate minimum in x-direction"
+{
+  : ::
+} -1.0
+REAL grid_ymin "Coordinate minimum in y-direction"
+{
+  : ::
+} -1.0
+REAL grid_zmin "Coordinate minimum in z-direction"
+{
+  : ::
+} -1.0
+REAL grid_xmax "Coordinate maximum in x-direction"
+{
+  : ::
+} 1.0
+REAL grid_ymax "Coordinate maximum in y-direction"
+{
+  : ::
+} 1.0
+REAL grid_zmax "Coordinate maximum in z-direction"
+{
+  : ::
+} 1.0
+
+
+KEYWORD grid "Grid type"
+{
+  "box" :: "Box"
+  "octant" :: "Octant"
+  "full" :: "Full"
+  "minmax" :: "Minmax"
+} "box"
+
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..964b37b
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,11 @@
+# Schedule definitions for thorn CartGrid3D
+# $Header$
+
+STORAGE: grid::coordinates
+
+schedule CartGrid3D at CACTUS_BASEGRID
+{
+  LANG:Fortran
+} "Set up spatial 3D Cartesian coordinates on the GH"
+
+
diff --git a/src/CartGrid3D.F b/src/CartGrid3D.F
new file mode 100644
index 0000000..5a9332d
--- /dev/null
+++ b/src/CartGrid3D.F
@@ -0,0 +1,173 @@
+#define THORN_IS_grid
+
+#include "cctk.h"
+#include "declare_arguments.h"
+#include "declare_parameters.h"
+
+      subroutine CartGrid3D(CCTK_FARGUMENTS)
+
+      implicit none
+
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_PARAMETERS
+
+      integer :: iflag, iconv, i, j, k
+      REAL :: x0,y0,z0, dx, dy, dz
+
+      integer :: CCTK_Equals
+
+      iflag = 0
+
+      iconv = 2**convlevel
+
+      print *,"grid_xmin = ",grid_xmin
+
+c     Defining all the parameters
+c     ---------------------------
+      grid_xmin = -1.0
+      grid_xmax = 1.0
+      grid_ymin = -1.0
+      grid_ymax = 1.0
+      grid_zmin = -1.0
+      grid_zmax = 1.0
+
+c     First, grids that ignore dx,dy,dz
+c     ---------------------------------
+
+c      if (CCTK_Equals(grid,'minmax')) then
+
+         iflag = iflag +1
+
+c        Set minimum values of coordinates
+         x0 = grid_xmin
+         y0 = grid_ymin
+         z0 = grid_zmin
+
+c        dx,dy,dz on the coarsest grid of each GH
+         coarse_dx = (grid_xmax-grid_xmin)/max(global_sh(1)-1,1)
+         coarse_dy = (grid_ymax-grid_ymin)/max(global_sh(2)-1,1)
+         coarse_dz = (grid_zmax-grid_zmin)/max(global_sh(3)-1,1)
+
+c        dx,dy,dz on the grid we are on
+         dx = coarse_dx/levfac
+         dy = coarse_dy/levfac
+         dz = coarse_dz/levfac
+
+         write(*,'(1X,A,1X,3(A,G12.7,2X))') 
+     &        'Minmax grid: ','dx=>',dx,'dy=>',dy,'dz=>',dz
+
+c      elseif (CCTK_Equals(grid,'box')) then 
+
+c         iflag = iflag + 1
+ 
+c         x0 = -0.5
+c         y0 = -0.5
+c         z0 = -0.5
+c         dx = 1.d0/max(global_nx-1,1)
+c         dy = 1.d0/max(global_ny-1,1)
+c         dz = 1.d0/max(global_nz-1,1)
+c         if (global_nx == 1) x0 = 0.0D0
+c         if (global_ny == 1) y0 = 0.0D0
+c         if (global_nz == 1) z0 = 0.0D0
+
+c         write(*,'(1X,A,1X,3(A,G12.7,2X))') 
+c     &        'Box grid: Resetting ','dx=>',dx,'dy=>',dy,'dz=>',dz
+
+c     Now, grids that use dx,dy,dz
+c     ----------------------------
+
+c      else
+
+c         dx = grid_dx*iconv
+c         dy = grid_dy*iconv
+c         dz = grid_dz*iconv
+
+c         if (CCTK_Equals(grid,'bitant')) then 
+
+c            iflag = iflag + 1
+c            x0 = (0.5 - global_nx/2)*dx
+c            y0 = (0.5 - global_ny/2)*dy
+
+c            if(nghostzones == 2) then
+c               z0 = -1.5*dz
+c            else
+c               z0 = -0.5*dz
+c            end if
+
+c         end if
+
+c         if (CCTK_Equals(grid,'quadrant')) then 
+
+c            iflag = iflag + 1
+
+c            if(nghostzones == 2) then
+c               x0 = -1.5D0*dx
+c               y0 = -1.5D0*dy
+c            else
+c               x0 = -0.5D0*dx
+c               y0 = -0.5D0*dy
+c            end if
+	    
+c            z0 = (0.5D0 - global_nz/2)*dz
+	    
+c	  end if
+
+c	  if (CCTK_Equals(grid,'octant')) then 
+	    
+c            if(nghostzones == 2) then
+	    
+c	      iflag = iflag + 1
+c	      x0 = -1.5*dx
+c	      y0 = -1.5*dy
+c	      z0 = -1.5*dz
+	    
+c            else
+	    
+c	      iflag = iflag + 1
+c	      x0 =  -0.5*dx
+c	      y0 =  -0.5*dy
+c	      z0 =  -0.5*dz
+	    
+c            end if
+
+c         end if
+
+c         if (CCTK_Equals(grid,'full')) then 
+
+c            iflag = iflag + 1
+c            x0 = (0.5 - global_nx/2)*dx
+c            y0 = (0.5 - global_ny/2)*dy
+c            z0 = (0.5 - global_nz/2)*dz
+
+c         end if
+
+c      end if
+
+
+c     No grid was set up
+c     ------------------
+c      if (iflag.ne.1) then
+c        write(*,*) 'setbasegrid error: check your grid keyword'
+c        STOP
+c     endif
+
+
+c     Set spatial coordinates
+c     -----------------------
+      do i=1,sh(1)
+         do j=1,sh(2)
+            do k=1,sh(3)
+               x(i,j,k) = dx*(i+lb(1)-1) + x0
+               y(i,j,k) = dy*(j+lb(2)-1) + y0
+               z(i,j,k) = dz*(k+lb(3)-1) + z0
+            end do
+         end do
+      end do
+
+      r = sqrt(x**2 + y**2 + z**2)
+
+      return
+      end subroutine CartGrid3D
+
+
+
diff --git a/src/WaveToy.F b/src/WaveToy.F
new file mode 100644
index 0000000..f23ec94
--- /dev/null
+++ b/src/WaveToy.F
@@ -0,0 +1,192 @@
+#include "cctk.h"
+#include "declare_parameters.h"
+#include "declare_arguments.h"
+
+      subroutine WaveToy_Boundaries(CCTK_FARGUMENTS)
+
+      implicit none
+
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_PARAMETERS
+
+      INTEGER CCTK_Equals
+      
+      WaveToy_bound = "none"
+
+      if(.NOT.CCTK_Equals(WaveToy_Bound, "periodic") .and. 
+     &     .NOT.CCTK_Equals(WaveToy_Bound,"none"))  then
+         scalar_field(1,:,:) = 0.0
+         scalar_field(:,1,:) = 0.0
+         scalar_field(:,:,1) = 0.0
+         
+         scalar_field(sh(1),:,:) = 0.0
+         scalar_field(:,sh(2),:) = 0.0
+         scalar_field(:,:,sh(3)) = 0.0
+      end if
+
+c     Should have a call to symmetry bounds here
+
+      end subroutine WaveToy_boundaries
+
+
+
+      subroutine WaveToy_evolution(CCTK_FARGUMENTS)
+
+      implicit none
+
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_PARAMETERS
+
+c     Declare local variables
+      INTEGER :: i,j,k
+      INTEGER :: istart, jstart, kstart, iend, jend, kend
+      REAL :: dx,dy,dz
+
+      dx = coarse_dx/levfac
+      dy = coarse_dy/levfac
+      dz = coarse_dz/levfac
+
+      istart = 2
+      jstart = 2
+      kstart = 2
+      
+      iend = sh(1)-1
+      jend = sh(2)-1
+      kend = sh(3)-1
+
+      do k = kstart, kend
+         do j = jstart, jend
+            do i = istart, iend
+               
+               scalar_field_tmp(i,j,k) = 
+     1		 2.0*(1.0 - (dt**2)*(1.0/dx**2 +   
+     2              1.0/dy**2 +1.0/dz**2))*scalar_field(i,j,k) -
+     3              scalar_field_old(i,j,k) +
+     4              (dt**2) * 
+     5              ((scalar_field(i+1,j,k)+scalar_field(i-1,j,k))/dx**2
+     6              +(scalar_field(i,j+1,k)+scalar_field(i,j-1,k))/dy**2
+     7              +(scalar_field(i,j,k+1)+scalar_field(i,j,k-1))/dz**2)
+
+            end do
+         end do
+      end do
+         
+      scalar_field_old = scalar_field
+      scalar_field     = scalar_field_tmp
+      
+      call WaveToy_boundaries(CCTK_FARGUMENTS)
+
+c      call CCTK_SyncGroup("scalarfields")
+
+      end subroutine WaveToy_evolution
+
+
+
+      subroutine WaveToy_init(CCTK_FARGUMENTS)
+
+      implicit none
+
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_PARAMETERS
+
+      REAL :: omega, pi, kx, ky, kz, radius, amplitude
+
+      integer CCTK_Equals
+
+      pi = 4.0*atan(1.0)
+
+c     Set parameters manually for now
+c     -------------------------------
+      WaveToy_dtfac = 0.5
+      WaveToy_kx = 2
+      WaveToy_ky = 2
+      WaveToy_kz = 2
+      WaveToy_radius = 0.25
+      WaveToy_amplitude = 1.0
+      WaveToy_initialdata = "plane"
+
+c     Calculate timestep
+      dt        = WaveToy_dtfac*coarse_dx/levfac
+      time = 0
+
+c     Shorthand
+      kx        = WaveToy_kx
+      ky        = WaveToy_ky
+      kz        = WaveToy_kz
+      radius    = WaveToy_radius
+      amplitude = WaveToy_amplitude
+
+      call WaveToy_boundaries(CCTK_FARGUMENTS)
+
+c      call CCTK_SyncGroup("scalarfields")
+
+      omega= sqrt(kx**2+ky**2+kz**2)
+
+c     Initialise the scalar field.  Can be plane waves, spherical waves
+c     or waves in a box.
+
+      if (CCTK_Equals(WaveToy_InitialData,"plane")) then 
+         scalar_field     = amplitude*
+     &        cos(kx*x+ky*y+kz*z+omega*time)
+         scalar_field_old = amplitude*
+     &        cos(kx*x+ky*y+kz*z+omega*(time-dt))
+
+      else if (CCTK_Equals(WaveToy_initialdata,"spherical")) then 
+
+         scalar_field     = amplitude*cos(kx*sqrt(x**2+y**2+z**2)+omega*time)
+         scalar_field_old = amplitude*cos(kx*sqrt(x**2+y**2+z**2)+omega*(time-dt))
+
+      else if (CCTK_Equals(WaveToy_InitialData, "box")) then
+
+c     Only currently implemented for the "box" grid type.
+
+c         if (CCTK_Equals(grid, "box")) then
+
+c      Must have all the wavenumbers non-zero for this to work.
+
+c            if (kx == 0) then
+c               write(*,*) "Box mode selected.  Setting kx to 1"
+c               kx = 1
+c            end if
+
+c            if (ky == 0) then
+c               write(*,*) "Box mode selected.  Setting ky to 1"
+c               ky = 1
+c            end if
+
+c            if (kz == 0) then
+c               write(*,*) "Box mode selected.  Setting kz to 1"
+c               kz = 1
+c            end if
+
+c     Set up the initial data.  
+c     Use kx,ky,kz as number of modes in each direction.
+
+c            scalar_field = amplitude*sin(kx*(x-0.5)*pi)*
+c     $           sin(ky*(y-0.5)*pi)*
+c     $           sin(kz*(z-0.5)*pi)*
+c     $           cos(omega*time*pi)
+            
+c            scalar_field_old = amplitude*sin(kx*(x-0.5)*pi)*
+c     $           sin(ky*(y-0.5)*pi)*
+c     $           sin(kz*(z-0.5)*pi)*
+c     $           cos(omega*(time-dt)*pi)
+c         else
+c            write(*,*) "WaveToy box mode selected, but grid not set to box"
+c            STOP
+c         end if
+
+      else if (CCTK_Equals(WaveToy_InitialData,"tsgaussian")) then
+
+         scalar_field = exp(-amplitude*(sqrt(x**2+y**2+z**2)-radius)**2)
+         scalar_field_old = scalar_field
+
+      else
+
+         print *,"Unknown initialisation mode in WaveToy"
+c         CCTK_Warning(GH,0,"Unknown initialisation mode in WaveToy")
+c         CCTK_Stop(GH,0)
+
+      end if
+
+      end subroutine WaveToy_Init
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..bf2ffff
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn CartGrid3D
+# $Header$
+
+# Source files in this directory
+SRCS = CartGrid3D.F
+
+# Subdirectories containing source files
+SUBDIRS = 
+