Initial revision

darcs-hash:20010301114010-f6438-12fb8a9ffcc80e86c0a97e37b5b0dae0dbc59b79.gz

12 files changed, 435 insertions, 0 deletions

diff --git a/CarpetExtra/IDHydroToy/README b/CarpetExtra/IDHydroToy/README
new file mode 100644
index 000000000..3ab11b194
--- /dev/null
+++ b/CarpetExtra/IDHydroToy/README
@@ -0,0 +1,7 @@
+Cactus Code Thorn IDHydroToy
+Authors    : Erik Schnetter <schnetter@uni-tuebingen.de>
+CVS info   : $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/README,v 1.1 2001/03/18 22:37:10 eschnett Exp $
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
diff --git a/CarpetExtra/IDHydroToy/interface.ccl b/CarpetExtra/IDHydroToy/interface.ccl
new file mode 100644
index 000000000..2e39f2268
--- /dev/null
+++ b/CarpetExtra/IDHydroToy/interface.ccl
@@ -0,0 +1,5 @@
+# Interface definition for thorn IDHydroToy
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/interface.ccl,v 1.1 2001/03/18 22:37:10 eschnett Exp $
+
+implements: idhydrotoy
+inherits: hydrotoy grid
diff --git a/CarpetExtra/IDHydroToy/param.ccl b/CarpetExtra/IDHydroToy/param.ccl
new file mode 100644
index 000000000..a0fd93f5b
--- /dev/null
+++ b/CarpetExtra/IDHydroToy/param.ccl
@@ -0,0 +1,46 @@
+# Parameter definitions for thorn IDHydroToy
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/param.ccl,v 1.2 2002/03/23 20:20:59 schnetter Exp $
+
+restricted:
+
+KEYWORD initial_data "Type of initial data"
+{
+  "plane"      :: "Plane wave"
+  "gaussian"   :: "Gaussian wave"
+  "box"        :: "Box wave"
+  "none"       :: "No initial data, zero phi"
+} "gaussian"
+
+private:
+
+## Parameter for initial wavepulses
+
+REAL radius "The radius of the gaussian wave"
+{
+ 0:*  :: "Positive"
+} 0.0
+
+REAL sigma "The sigma for the gaussian wave"
+{
+ 0:* :: "Positive"
+} 0.1
+
+REAL kx "The wave number in the x-direction"
+{
+ *:* :: "No restriction"
+} 4.0
+
+REAL ky "The wave number in the y-direction"
+{
+ *:* :: "No restriction"
+} 0.0
+
+REAL kz "The wave number in the z-direction"
+{
+ *:* :: "No restriction"
+} 0.0
+
+REAL amplitude "The amplitude of the waves"
+{
+  *:* :: "No restriction"
+} 1.0
diff --git a/CarpetExtra/IDHydroToy/schedule.ccl b/CarpetExtra/IDHydroToy/schedule.ccl
new file mode 100644
index 000000000..c47184aab
--- /dev/null
+++ b/CarpetExtra/IDHydroToy/schedule.ccl
@@ -0,0 +1,13 @@
+# Schedule definitions for thorn IDHydroToy
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/schedule.ccl,v 1.3 2003/11/05 16:18:40 schnetter Exp $
+
+schedule IDHydroToy_Startup at STARTUP
+{
+  LANG: Fortran
+} "Register banner"
+
+schedule IDHydroToy_InitialData as HydroToy_InitialData at INITIAL
+{
+  LANG: Fortran
+  STORAGE: hydrotoy::hydroevolve[3]
+} "Initial data for 3D wave equation"
diff --git a/CarpetExtra/IDHydroToy/src/InitialData.F77 b/CarpetExtra/IDHydroToy/src/InitialData.F77
new file mode 100644
index 000000000..65f6e4e3a
--- /dev/null
+++ b/CarpetExtra/IDHydroToy/src/InitialData.F77
@@ -0,0 +1,231 @@
+c     -*-Fortran-*-
+c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/src/InitialData.F77,v 1.6 2003/11/05 16:18:40 schnetter Exp $
+      
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+      
+      subroutine IDHydroToy_InitialData (CCTK_ARGUMENTS)
+      
+      implicit none
+      
+      DECLARE_CCTK_ARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+      
+      CCTK_REAL pi
+      CCTK_REAL omega
+      CCTK_REAL dt
+      CCTK_REAL x,y,z, r
+      integer i,j,k
+      
+      CCTK_REAL vr
+      
+      external erf
+      real*8 erf
+      
+      pi = 4*atan(1.d0)
+      
+      omega = sqrt(kx**2+ky**2+kz**2)
+      
+      dt = CCTK_DELTA_TIME
+      
+      if (CCTK_EQUALS(initial_data,"plane")) then 
+         
+         do k=1,cctk_lsh(3)
+            do j=1,cctk_lsh(2)
+               do i=1,cctk_lsh(1)
+                  
+                  x = cart3d_x(i,j,k)
+                  y = cart3d_y(i,j,k)
+                  z = cart3d_z(i,j,k)
+                  
+                  u(i,j,k) = amplitude
+     $                 * cos((kx*x + ky*y + kz*z + omega*cctk_time) * pi)
+                  vx(i,j,k) = u(i,j,k) * kx / omega
+                  vy(i,j,k) = u(i,j,k) * ky / omega
+                  vz(i,j,k) = u(i,j,k) * kz / omega
+                  
+                  u_p(i,j,k) = amplitude
+     $                 * cos((kx*x + ky*y + kz*z + omega*(cctk_time-dt)) * pi)
+                  vx_p(i,j,k) = u_p(i,j,k) * kx / omega
+                  vy_p(i,j,k) = u_p(i,j,k) * ky / omega
+                  vz_p(i,j,k) = u_p(i,j,k) * kz / omega
+                  
+                  u_p_p(i,j,k) = amplitude
+     $                 * cos((kx*x + ky*y + kz*z + omega*(cctk_time-2*dt)) * pi)
+                  vx_p_p(i,j,k) = u_p_p(i,j,k) * kx / omega
+                  vy_p_p(i,j,k) = u_p_p(i,j,k) * ky / omega
+                  vz_p_p(i,j,k) = u_p_p(i,j,k) * kz / omega
+                  
+               end do
+            end do
+         end do
+         
+      else if (CCTK_EQUALS(initial_data,"gaussian")) then 
+         
+         do k=1,cctk_lsh(3)
+            do j=1,cctk_lsh(2)
+               do i=1,cctk_lsh(1)
+                  
+                  x = cart3d_x(i,j,k)
+                  y = cart3d_y(i,j,k)
+                  z = cart3d_z(i,j,k)
+                  r = spher3d_r(i,j,k)
+                  
+                  u(i,j,k) = amplitude
+     $                 * exp(- (r - radius)**2 / sigma**2)
+                  
+                  vr = - 2*amplitude * (r - radius) / sigma**2
+     $                 * exp(- (r - radius)**2 / sigma**2)
+                  vx(i,j,k) = vr * x/r
+                  vy(i,j,k) = vr * y/r
+                  vz(i,j,k) = vr * z/r
+                  
+                  u_p(i,j,k) = amplitude/2 * (r - dt) / r
+     $                 * exp(- (r - radius - dt)**2 / sigma**2)
+     $                 + amplitude/2 * (r + dt) / r
+     $                 * exp(- (r - radius + dt)**2 / sigma**2)
+                  
+                  vr = - amplitude/2 * (-dt / r**2 + (r - dt) * (r - radius - dt) / (r * sigma**2))
+     $                 * exp(- (r - radius - dt)**2 / sigma**2)
+     $                 - amplitude/2 * ( dt / r**2 + (r + dt) * (r - radius + dt) / (r * sigma**2))
+     $                 * exp(- (r - radius - dt)**2 / sigma**2)
+                  vx_p(i,j,k) = vr * x/r
+                  vy_p(i,j,k) = vr * y/r
+                  vz_p(i,j,k) = vr * z/r
+                  
+                  u_p_p(i,j,k) = amplitude/2 * (r - 2*dt) / r
+     $                 * exp(- (r - radius - 2*dt)**2 / sigma**2)
+     $                 + amplitude/2 * (r + 2*dt) / r
+     $                 * exp(- (r - radius + 2*dt)**2 / sigma**2)
+                  
+                  vr = - amplitude/2 * (-2*dt / r**2 + (r - 2*dt) * (r - radius - 2*dt) / (r * sigma**2))
+     $                 * exp(- (r - radius - 2*dt)**2 / sigma**2)
+     $                 - amplitude/2 * ( 2*dt / r**2 + (r + 2*dt) * (r - radius + 2*dt) / (r * sigma**2))
+     $                 * exp(- (r - radius - 2*dt)**2 / sigma**2)
+                  vx_p_p(i,j,k) = vr * x/r
+                  vy_p_p(i,j,k) = vr * y/r
+                  vz_p_p(i,j,k) = vr * z/r
+                  
+               end do
+            end do
+         end do
+         
+      else if (CCTK_EQUALS(initial_data, "box")) then
+         
+c     Use kx,ky,kz as number of modes in each direction.
+         
+         do k=1,cctk_lsh(3)
+            do j=1,cctk_lsh(2)
+               do i=1,cctk_lsh(1)
+                  
+                  x = cart3d_x(i,j,k)
+                  y = cart3d_y(i,j,k)
+                  z = cart3d_z(i,j,k)
+                  
+                  u(i,j,k) = amplitude
+     $                 * sin(kx * (x - 0.5d0) * pi)
+     $                 * sin(ky * (y - 0.5d0) * pi)
+     $                 * sin(kz * (z - 0.5d0) * pi)
+     $                 * cos(omega * cctk_time * pi)
+                  vx(i,j,k) = amplitude
+     $                 * cos(kx * (x - 0.5d0) * pi)
+     $                 * sin(ky * (y - 0.5d0) * pi)
+     $                 * sin(kz * (z - 0.5d0) * pi)
+     $                 * sin(omega * cctk_time * pi)
+     $                 * kx / omega
+                  vy(i,j,k) = amplitude
+     $                 * sin(kx * (x - 0.5d0) * pi)
+     $                 * cos(ky * (y - 0.5d0) * pi)
+     $                 * sin(kz * (z - 0.5d0) * pi)
+     $                 * sin(omega * cctk_time * pi)
+     $                 * ky / omega
+                  vz(i,j,k) = amplitude
+     $                 * sin(kx * (x - 0.5d0) * pi)
+     $                 * sin(ky * (y - 0.5d0) * pi)
+     $                 * cos(kz * (z - 0.5d0) * pi)
+     $                 * sin(omega * cctk_time * pi)
+     $                 * kz / omega
+                  
+                  u_p(i,j,k) = amplitude
+     $                 * sin(kx * (x - 0.5d0) * pi)
+     $                 * sin(ky * (y - 0.5d0) * pi)
+     $                 * sin(kz * (z - 0.5d0) * pi)
+     $                 * cos(omega * (cctk_time - dt) * pi)
+                  vx_p(i,j,k) = amplitude
+     $                 * cos(kx * (x - 0.5d0) * pi)
+     $                 * sin(ky * (y - 0.5d0) * pi)
+     $                 * sin(kz * (z - 0.5d0) * pi)
+     $                 * sin(omega * (cctk_time - dt) * pi)
+     $                 * kx / omega
+                  vy_p(i,j,k) = amplitude
+     $                 * sin(kx * (x - 0.5d0) * pi)
+     $                 * cos(ky * (y - 0.5d0) * pi)
+     $                 * sin(kz * (z - 0.5d0) * pi)
+     $                 * sin(omega * (cctk_time - dt) * pi)
+     $                 * ky / omega
+                  vz_p(i,j,k) = amplitude
+     $                 * sin(kx * (x - 0.5d0) * pi)
+     $                 * sin(ky * (y - 0.5d0) * pi)
+     $                 * cos(kz * (z - 0.5d0) * pi)
+     $                 * sin(omega * (cctk_time - dt) * pi)
+     $                 * kz / omega
+                  
+                  u_p_p(i,j,k) = amplitude
+     $                 * sin(kx * (x - 0.5d0) * pi)
+     $                 * sin(ky * (y - 0.5d0) * pi)
+     $                 * sin(kz * (z - 0.5d0) * pi)
+     $                 * cos(omega * (cctk_time - 2*dt) * pi)
+                  vx_p_p(i,j,k) = amplitude
+     $                 * cos(kx * (x - 0.5d0) * pi)
+     $                 * sin(ky * (y - 0.5d0) * pi)
+     $                 * sin(kz * (z - 0.5d0) * pi)
+     $                 * sin(omega * (cctk_time - 2*dt) * pi)
+     $                 * kx / omega
+                  vy_p_p(i,j,k) = amplitude
+     $                 * sin(kx * (x - 0.5d0) * pi)
+     $                 * cos(ky * (y - 0.5d0) * pi)
+     $                 * sin(kz * (z - 0.5d0) * pi)
+     $                 * sin(omega * (cctk_time - 2*dt) * pi)
+     $                 * ky / omega
+                  vz_p_p(i,j,k) = amplitude
+     $                 * sin(kx * (x - 0.5d0) * pi)
+     $                 * sin(ky * (y - 0.5d0) * pi)
+     $                 * cos(kz * (z - 0.5d0) * pi)
+     $                 * sin(omega * (cctk_time - 2*dt) * pi)
+     $                 * kz / omega
+                  
+               end do
+            end do
+         end do
+         
+      else
+         
+         do k=1,cctk_lsh(3)
+            do j=1,cctk_lsh(2)
+               do i=1,cctk_lsh(1)
+                  
+                  u(i,j,k)  = 0
+                  vx(i,j,k) = 0
+                  vy(i,j,k) = 0
+                  vz(i,j,k) = 0
+                  
+                  u_p(i,j,k)  = 0
+                  vx_p(i,j,k) = 0
+                  vy_p(i,j,k) = 0
+                  vz_p(i,j,k) = 0
+                  
+                  u_p_p(i,j,k)  = 0
+                  vx_p_p(i,j,k) = 0
+                  vy_p_p(i,j,k) = 0
+                  vz_p_p(i,j,k) = 0
+                  
+               end do
+            end do
+         end do 
+         
+      end if
+      
+      end 
diff --git a/CarpetExtra/IDHydroToy/src/Startup.F77 b/CarpetExtra/IDHydroToy/src/Startup.F77
new file mode 100644
index 000000000..d4e87f0af
--- /dev/null
+++ b/CarpetExtra/IDHydroToy/src/Startup.F77
@@ -0,0 +1,15 @@
+c     -*-Fortran-*-
+c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/src/Startup.F77,v 1.1 2001/03/18 22:37:10 eschnett Exp $
+      
+#include "cctk.h"
+      
+      subroutine IDHydroToy_Startup
+      
+      implicit none
+      
+      integer ierr
+      
+      call CCTK_RegisterBanner
+     $     (ierr, "IDHydroToy: Initial data for a Scalar and a Vector Field")
+      
+      end 
diff --git a/CarpetExtra/IDHydroToy/src/erf.f77 b/CarpetExtra/IDHydroToy/src/erf.f77
new file mode 100644
index 000000000..30bc84520
--- /dev/null
+++ b/CarpetExtra/IDHydroToy/src/erf.f77
@@ -0,0 +1,13 @@
+      FUNCTION erf(x)
+      implicit none
+      REAL*8 erf,x
+CU    USES gammp
+      REAL*8 gammp
+      if(x.lt.0)then
+        erf=-gammp(.5d0,x**2)
+      else
+        erf=gammp(.5d0,x**2)
+      endif
+      return
+      END
+C  (C) Copr. 1986-92 Numerical Recipes Software t4.
diff --git a/CarpetExtra/IDHydroToy/src/gammln.f77 b/CarpetExtra/IDHydroToy/src/gammln.f77
new file mode 100644
index 000000000..b99fb6aa2
--- /dev/null
+++ b/CarpetExtra/IDHydroToy/src/gammln.f77
@@ -0,0 +1,22 @@
+      FUNCTION gammln(xx)
+      implicit none
+      REAL*8 gammln,xx
+      INTEGER j
+      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
+      SAVE cof,stp
+      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
+     *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
+     *-.5395239384953d-5,2.5066282746310005d0/
+      x=xx
+      y=x
+      tmp=x+5.5d0
+      tmp=(x+0.5d0)*log(tmp)-tmp
+      ser=1.000000000190015d0
+      do 11 j=1,6
+        y=y+1.d0
+        ser=ser+cof(j)/y
+11    continue
+      gammln=tmp+log(stp*ser/x)
+      return
+      END
+C  (C) Copr. 1986-92 Numerical Recipes Software t4.
diff --git a/CarpetExtra/IDHydroToy/src/gammp.f77 b/CarpetExtra/IDHydroToy/src/gammp.f77
new file mode 100644
index 000000000..99a938857
--- /dev/null
+++ b/CarpetExtra/IDHydroToy/src/gammp.f77
@@ -0,0 +1,16 @@
+      FUNCTION gammp(a,x)
+      implicit none
+      REAL*8 a,gammp,x
+CU    USES gcf,gser
+      REAL*8 gammcf,gamser,gln
+      if(x.lt.0.or.a.le.0)pause 'bad arguments in gammp'
+      if(x.lt.a+1)then
+        call gser(gamser,a,x,gln)
+        gammp=gamser
+      else
+        call gcf(gammcf,a,x,gln)
+        gammp=1-gammcf
+      endif
+      return
+      END
+C  (C) Copr. 1986-92 Numerical Recipes Software t4.
diff --git a/CarpetExtra/IDHydroToy/src/gcf.f77 b/CarpetExtra/IDHydroToy/src/gcf.f77
new file mode 100644
index 000000000..09fd58aa0
--- /dev/null
+++ b/CarpetExtra/IDHydroToy/src/gcf.f77
@@ -0,0 +1,30 @@
+      SUBROUTINE gcf(gammcf,a,x,gln)
+      implicit none
+      INTEGER ITMAX
+      REAL*8 a,gammcf,gln,x,EPS,FPMIN
+      PARAMETER (ITMAX=100,EPS=3.d-7,FPMIN=1.d-30)
+CU    USES gammln
+      INTEGER i
+      REAL*8 an,b,c,d,del,h,gammln
+      gln=gammln(a)
+      b=x+1-a
+      c=1/FPMIN
+      d=1/b
+      h=d
+      do 11 i=1,ITMAX
+        an=-i*(i-a)
+        b=b+2
+        d=an*d+b
+        if(abs(d).lt.FPMIN)d=FPMIN
+        c=b+an/c
+        if(abs(c).lt.FPMIN)c=FPMIN
+        d=1/d
+        del=d*c
+        h=h*del
+        if(abs(del-1).lt.EPS)goto 1
+11    continue
+      pause 'a too large, ITMAX too small in gcf'
+1     gammcf=exp(-x+a*log(x)-gln)*h
+      return
+      END
+C  (C) Copr. 1986-92 Numerical Recipes Software t4.
diff --git a/CarpetExtra/IDHydroToy/src/gser.f77 b/CarpetExtra/IDHydroToy/src/gser.f77
new file mode 100644
index 000000000..fad2be01f
--- /dev/null
+++ b/CarpetExtra/IDHydroToy/src/gser.f77
@@ -0,0 +1,28 @@
+      SUBROUTINE gser(gamser,a,x,gln)
+      implicit none
+      INTEGER ITMAX
+      REAL*8 a,gamser,gln,x,EPS
+      PARAMETER (ITMAX=100,EPS=3.d-7)
+CU    USES gammln
+      INTEGER n
+      REAL*8 ap,del,sum,gammln
+      gln=gammln(a)
+      if(x.le.0)then
+        if(x.lt.0)pause 'x < 0 in gser'
+        gamser=0
+        return
+      endif
+      ap=a
+      sum=1/a
+      del=sum
+      do 11 n=1,ITMAX
+        ap=ap+1
+        del=del*x/ap
+        sum=sum+del
+        if(abs(del).lt.abs(sum)*EPS)goto 1
+11    continue
+      pause 'a too large, ITMAX too small in gser'
+1     gamser=sum*exp(-x+a*log(x)-gln)
+      return
+      END
+C  (C) Copr. 1986-92 Numerical Recipes Software t4.
diff --git a/CarpetExtra/IDHydroToy/src/make.code.defn b/CarpetExtra/IDHydroToy/src/make.code.defn
new file mode 100644
index 000000000..a63bda2f5
--- /dev/null
+++ b/CarpetExtra/IDHydroToy/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn IDHydroToy
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/src/make.code.defn,v 1.3 2003/09/20 13:47:27 schnetter Exp $
+
+# Source files in this directory
+SRCS = InitialData.F77 Startup.F77 erf.f77 gammln.f77 gammp.f77 gcf.f77 gser.f77
+
+# Subdirectories containing source files
+SUBDIRS = 
+