aboutsummaryrefslogtreecommitdiff
path: root/CarpetExtra/IDHydroToy/src/InitialData.F77
diff options
context:
space:
mode:
Diffstat (limited to 'CarpetExtra/IDHydroToy/src/InitialData.F77')
-rw-r--r--CarpetExtra/IDHydroToy/src/InitialData.F77163
1 files changed, 22 insertions, 141 deletions
diff --git a/CarpetExtra/IDHydroToy/src/InitialData.F77 b/CarpetExtra/IDHydroToy/src/InitialData.F77
index 65f6e4e3a..7c03722ad 100644
--- a/CarpetExtra/IDHydroToy/src/InitialData.F77
+++ b/CarpetExtra/IDHydroToy/src/InitialData.F77
@@ -1,10 +1,9 @@
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 $
+c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/src/InitialData.F77,v 1.1 2001/03/18 22:37:10 eschnett Exp $
#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Functions.h"
#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
subroutine IDHydroToy_InitialData (CCTK_ARGUMENTS)
@@ -16,49 +15,25 @@ c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/src/Initial
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)
+ $ * cos((kx*cart3d_x(i,j,k) + ky*cart3d_y(i,j,k)
+ $ + kz*cart3d_z(i,j,k) + 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
@@ -69,45 +44,12 @@ c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/src/Initial
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
+ $ * exp(- (spher3d_r(i,j,k) - radius + cctk_time)**2
+ $ / sigma**2)
+ vx(i,j,k) = u(i,j,k) * cart3d_x(i,j,k)/spher3d_r(i,j,k)
+ vy(i,j,k) = u(i,j,k) * cart3d_y(i,j,k)/spher3d_r(i,j,k)
+ vz(i,j,k) = u(i,j,k) * cart3d_z(i,j,k)/spher3d_r(i,j,k)
end do
end do
@@ -121,82 +63,31 @@ c Use kx,ky,kz as number of modes in each direction.
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)
+ $ * sin(kx * (cart3d_x(i,j,k) - 0.5d0) * pi)
+ $ * sin(ky * (cart3d_y(i,j,k) - 0.5d0) * pi)
+ $ * sin(kz * (cart3d_z(i,j,k) - 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)
+ $ * cos(kx * (cart3d_x(i,j,k) - 0.5d0) * pi)
+ $ * sin(ky * (cart3d_y(i,j,k) - 0.5d0) * pi)
+ $ * sin(kz * (cart3d_z(i,j,k) - 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(kx * (cart3d_x(i,j,k) - 0.5d0) * pi)
+ $ * cos(ky * (cart3d_y(i,j,k) - 0.5d0) * pi)
+ $ * sin(kz * (cart3d_z(i,j,k) - 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(kx * (cart3d_x(i,j,k) - 0.5d0) * pi)
+ $ * sin(ky * (cart3d_y(i,j,k) - 0.5d0) * pi)
+ $ * cos(kz * (cart3d_z(i,j,k) - 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
@@ -212,16 +103,6 @@ c Use kx,ky,kz as number of modes in each direction.
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
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-02 23:17:54 +0000