From 310f0ea48d18866b773136aed11200b6eda6378b Mon Sep 17 00:00:00 2001 From: eschnett <> Date: Thu, 1 Mar 2001 11:40:00 +0000 Subject: Initial revision darcs-hash:20010301114010-f6438-12fb8a9ffcc80e86c0a97e37b5b0dae0dbc59b79.gz --- CarpetExtra/IDHydroToy/README | 7 + CarpetExtra/IDHydroToy/interface.ccl | 5 + CarpetExtra/IDHydroToy/param.ccl | 46 ++++++ CarpetExtra/IDHydroToy/schedule.ccl | 13 ++ CarpetExtra/IDHydroToy/src/InitialData.F77 | 231 +++++++++++++++++++++++++++++ CarpetExtra/IDHydroToy/src/Startup.F77 | 15 ++ CarpetExtra/IDHydroToy/src/erf.f77 | 13 ++ CarpetExtra/IDHydroToy/src/gammln.f77 | 22 +++ CarpetExtra/IDHydroToy/src/gammp.f77 | 16 ++ CarpetExtra/IDHydroToy/src/gcf.f77 | 30 ++++ CarpetExtra/IDHydroToy/src/gser.f77 | 28 ++++ CarpetExtra/IDHydroToy/src/make.code.defn | 9 ++ 12 files changed, 435 insertions(+) create mode 100644 CarpetExtra/IDHydroToy/README create mode 100644 CarpetExtra/IDHydroToy/interface.ccl create mode 100644 CarpetExtra/IDHydroToy/param.ccl create mode 100644 CarpetExtra/IDHydroToy/schedule.ccl create mode 100644 CarpetExtra/IDHydroToy/src/InitialData.F77 create mode 100644 CarpetExtra/IDHydroToy/src/Startup.F77 create mode 100644 CarpetExtra/IDHydroToy/src/erf.f77 create mode 100644 CarpetExtra/IDHydroToy/src/gammln.f77 create mode 100644 CarpetExtra/IDHydroToy/src/gammp.f77 create mode 100644 CarpetExtra/IDHydroToy/src/gcf.f77 create mode 100644 CarpetExtra/IDHydroToy/src/gser.f77 create mode 100644 CarpetExtra/IDHydroToy/src/make.code.defn (limited to 'CarpetExtra/IDHydroToy') 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 +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 = + -- cgit v1.2.3