Using three time levels.

darcs-hash:20010326002850-f6438-82257c36db375047ccc8fa1ea38cf6e0a18ed30c.gz

3 files changed, 161 insertions, 24 deletions

diff --git a/CarpetExtra/IDHydroToy/src/InitialData.F77 b/CarpetExtra/IDHydroToy/src/InitialData.F77
index 5bbb98032..92bcf0e2a 100644
--- a/CarpetExtra/IDHydroToy/src/InitialData.F77
+++ b/CarpetExtra/IDHydroToy/src/InitialData.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/src/InitialData.F77,v 1.2 2001/03/19 21:30:37 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/src/InitialData.F77,v 1.3 2001/03/26 02:28:50 eschnett Exp $
       
 #include "cctk.h"
 #include "cctk_Parameters.h"
@@ -15,6 +15,7 @@ c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/src/Initial
       
       CCTK_REAL pi
       CCTK_REAL omega
+      CCTK_REAL dt
       integer i,j,k
       
       CCTK_REAL vr
@@ -26,6 +27,8 @@ c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/src/Initial
       
       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)
@@ -39,6 +42,13 @@ c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/src/Initial
                   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*cart3d_x(i,j,k) + ky*cart3d_y(i,j,k)
+     $                 + kz*cart3d_z(i,j,k) + 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
+                  
                end do
             end do
          end do
@@ -51,8 +61,7 @@ c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDHydroToy/src/Initial
                   
                   u(i,j,k) = amplitude / spher3d_r(i,j,k)
      $                 * exp(- (spher3d_r(i,j,k) - radius + cctk_time)**2 / sigma**2)
-
-c     vr = -1/2*A*(sqrt(Pi)*sigma*erf((t-r+R)/sigma)+2*exp(-(t-r+R)^2/(sigma^2))*r)/(r^2)
+                  
                   vr = -0.5d0 * amplitude
      $                 * (sqrt(pi) * sigma * erf((cctk_time - spher3d_r(i,j,k) + radius) / sigma)
      $                    + 2*exp(-(cctk_time - spher3d_r(i,j,k) + radius)**2 / sigma**2) * spher3d_r(i,j,k))
@@ -61,6 +70,17 @@ c     vr = -1/2*A*(sqrt(Pi)*sigma*erf((t-r+R)/sigma)+2*exp(-(t-r+R)^2/(sigma^2))
                   vy(i,j,k) = vr * cart3d_y(i,j,k)/spher3d_r(i,j,k)
                   vz(i,j,k) = vr * cart3d_z(i,j,k)/spher3d_r(i,j,k)
                   
+                  u_p(i,j,k) = amplitude / spher3d_r(i,j,k)
+     $                 * exp(- (spher3d_r(i,j,k) - radius + (cctk_time-dt))**2 / sigma**2)
+                  
+                  vr = -0.5d0 * amplitude
+     $                 * (sqrt(pi) * sigma * erf((cctk_time - spher3d_r(i,j,k) + radius) / sigma)
+     $                    + 2*exp(-((cctk_time-dt) - spher3d_r(i,j,k) + radius)**2 / sigma**2) * spher3d_r(i,j,k))
+     $                 / spher3d_r(i,j,k)**2
+                  vx_p(i,j,k) = vr * cart3d_x(i,j,k)/spher3d_r(i,j,k)
+                  vy_p(i,j,k) = vr * cart3d_y(i,j,k)/spher3d_r(i,j,k)
+                  vz_p(i,j,k) = vr * cart3d_z(i,j,k)/spher3d_r(i,j,k)
+                  
                end do
             end do
          end do
@@ -78,7 +98,6 @@ c     Use kx,ky,kz as number of modes in each direction.
      $                 * 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 * (cart3d_x(i,j,k) - 0.5d0) * pi)
      $                 * sin(ky * (cart3d_y(i,j,k) - 0.5d0) * pi)
@@ -98,6 +117,30 @@ c     Use kx,ky,kz as number of modes in each direction.
      $                 * sin(omega * cctk_time * pi)
      $                 * kz / omega
                   
+                  u_p(i,j,k) = amplitude
+     $                 * 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-dt) * pi)
+                  vx_p(i,j,k) = amplitude
+     $                 * 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-dt) * pi)
+     $                 * kx / omega
+                  vy_p(i,j,k) = amplitude
+     $                 * 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-dt) * pi)
+     $                 * ky / omega
+                  vz_p(i,j,k) = amplitude
+     $                 * 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-dt) * pi)
+     $                 * kz / omega
+                  
                end do
             end do
          end do
@@ -113,6 +156,11 @@ 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
+                  
                end do
             end do
          end do 
diff --git a/CarpetExtra/IDSpaceTimeToy/src/InitialData.F77 b/CarpetExtra/IDSpaceTimeToy/src/InitialData.F77
index 4fcc252de..5582ed669 100644
--- a/CarpetExtra/IDSpaceTimeToy/src/InitialData.F77
+++ b/CarpetExtra/IDSpaceTimeToy/src/InitialData.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDSpaceTimeToy/src/InitialData.F77,v 1.3 2001/03/20 17:27:26 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDSpaceTimeToy/src/InitialData.F77,v 1.4 2001/03/26 02:28:54 eschnett Exp $
       
 #include "cctk.h"
 #include "cctk_Parameters.h"
@@ -15,12 +15,15 @@ c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDSpaceTimeToy/src/Ini
       
       CCTK_REAL pi
       CCTK_REAL omega
+      CCTK_REAL dt
       integer i,j,k
       
       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)
@@ -35,6 +38,14 @@ c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDSpaceTimeToy/src/Ini
      $                 + kz*cart3d_z(i,j,k) + omega*cctk_time) * pi)
      $                 * pi * omega
                   
+                  phi_p(i,j,k) = amplitude
+     $                 * cos((kx*cart3d_x(i,j,k) + ky*cart3d_y(i,j,k)
+     $                 + kz*cart3d_z(i,j,k) + omega*(cctk_time-dt)) * pi)
+                  psi_p(i,j,k) = - amplitude
+     $                 * sin((kx*cart3d_x(i,j,k) + ky*cart3d_y(i,j,k)
+     $                 + kz*cart3d_z(i,j,k) + omega*(cctk_time-dt)) * pi)
+     $                 * pi * omega
+                  
                end do
             end do
          end do
@@ -50,6 +61,11 @@ c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDSpaceTimeToy/src/Ini
                   psi(i,j,k) = phi(i,j,k)
      $                 * 2 * (spher3d_r(i,j,k) - radius - cctk_time) / sigma**2
                   
+                  phi_p(i,j,k) = amplitude / spher3d_r(i,j,k)
+     $                 * exp(- (spher3d_r(i,j,k) - radius - (cctk_time-dt))**2 / sigma**2)
+                  psi_p(i,j,k) = phi(i,j,k)
+     $                 * 2 * (spher3d_r(i,j,k) - radius - (cctk_time-dt)) / sigma**2
+                  
                end do
             end do
          end do
@@ -67,7 +83,6 @@ c     Use kx,ky,kz as number of modes in each direction.
      $                 * 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)
-                  
                   psi(i,j,k) = - amplitude
      $                 * sin(kx * (cart3d_x(i,j,k) - 0.5d0) * pi)
      $                 * sin(ky * (cart3d_y(i,j,k) - 0.5d0) * pi)
@@ -75,6 +90,18 @@ c     Use kx,ky,kz as number of modes in each direction.
      $                 * sin(omega * cctk_time * pi)
      $                 * omega * pi
                   
+                  phi_p(i,j,k) = amplitude
+     $                 * 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-dt) * pi)
+                  psi_p(i,j,k) = - amplitude
+     $                 * 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)
+     $                 * sin(omega * (cctk_time-dt) * pi)
+     $                 * omega * pi
+                  
                end do
             end do
          end do
@@ -88,20 +115,28 @@ c     Use kx,ky,kz as number of modes in each direction.
                   phi(i,j,k) = 0
                   psi(i,j,k) = 0
                   
+                  phi_p(i,j,k) = 0
+                  psi_p(i,j,k) = 0
+                  
                end do
             end do
          end do 
          
       end if
       
-      do k=1,cctk_lsh(3)
-         do j=1,cctk_lsh(2)
-            do i=1,cctk_lsh(1)
-               
-               psi(i,j,k) = psi(i,j,k) - u(i,j,k)
-               
+      if (cctk_isthornactive("hydrotoy").eq.1) then
+         
+         do k=1,cctk_lsh(3)
+            do j=1,cctk_lsh(2)
+               do i=1,cctk_lsh(1)
+                  
+                  psi(i,j,k) = psi(i,j,k) - u(i,j,k)
+                  psi_p(i,j,k) = psi_p(i,j,k) - u_p(i,j,k)
+                  
+               end do
             end do
          end do
-      end do 
+         
+      end if
       
       end 
diff --git a/CarpetExtra/SpaceTimeToy/src/SpaceTimeToy.F77 b/CarpetExtra/SpaceTimeToy/src/SpaceTimeToy.F77
index b9fb12aaf..ecb245028 100644
--- a/CarpetExtra/SpaceTimeToy/src/SpaceTimeToy.F77
+++ b/CarpetExtra/SpaceTimeToy/src/SpaceTimeToy.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/SpaceTimeToy/src/SpaceTimeToy.F77,v 1.6 2001/03/22 18:42:06 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/SpaceTimeToy/src/SpaceTimeToy.F77,v 1.7 2001/03/26 02:29:24 eschnett Exp $
       
 #include "cctk.h" 
 #include "cctk_Parameters.h"
@@ -25,15 +25,42 @@ c     Copy
                phi_i(i,j,k) = phi(i,j,k)
                psi_i(i,j,k) = psi(i,j,k)
                
-               u_i(i,j,k)  = u(i,j,k)
-               vx_i(i,j,k) = vx(i,j,k)
-               vy_i(i,j,k) = vy(i,j,k)
-               vz_i(i,j,k) = vz(i,j,k)
-               
             end do
          end do
       end do
       
+      if (cctk_isthornactive("hydrotoy").eq.1) then
+         
+         do k=1,cctk_lsh(3)
+            do j=1,cctk_lsh(2)
+               do i=1,cctk_lsh(1)
+               
+                  u_i(i,j,k)  = u(i,j,k)
+                  vx_i(i,j,k) = vx(i,j,k)
+                  vy_i(i,j,k) = vy(i,j,k)
+                  vz_i(i,j,k) = vz(i,j,k)
+                  
+               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(i,j,k)  = 0
+                  vx_i(i,j,k) = 0
+                  vy_i(i,j,k) = 0
+                  vz_i(i,j,k) = 0
+                  
+               end do
+            end do
+         end do
+         
+      end if
+      
 c     Evolve
       call SpaceTimeToy_Step (CCTK_PASS_FTOF)
       
@@ -72,15 +99,42 @@ c     Average
                phi_i(i,j,k) = half * (phi(i,j,k) + phi_n(i,j,k))
                psi_i(i,j,k) = half * (psi(i,j,k) + psi_n(i,j,k))
                
-               u_i(i,j,k)  = half * (u(i,j,k)  + u_n(i,j,k))
-               vx_i(i,j,k) = half * (vx(i,j,k) + vx_n(i,j,k))
-               vy_i(i,j,k) = half * (vy(i,j,k) + vy_n(i,j,k))
-               vz_i(i,j,k) = half * (vz(i,j,k) + vz_n(i,j,k))
-               
             end do
          end do
       end do
       
+      if (cctk_isthornactive("hydrotoy").eq.1) then
+         
+         do k=1,cctk_lsh(3)
+            do j=1,cctk_lsh(2)
+               do i=1,cctk_lsh(1)
+                  
+                  u_i(i,j,k)  = half * (u(i,j,k)  + u_n(i,j,k))
+                  vx_i(i,j,k) = half * (vx(i,j,k) + vx_n(i,j,k))
+                  vy_i(i,j,k) = half * (vy(i,j,k) + vy_n(i,j,k))
+                  vz_i(i,j,k) = half * (vz(i,j,k) + vz_n(i,j,k))
+                  
+               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(i,j,k)  = 0
+                  vx_i(i,j,k) = 0
+                  vy_i(i,j,k) = 0
+                  vz_i(i,j,k) = 0
+                  
+               end do
+            end do
+         end do
+         
+      end if
+      
 c     Evolve
       call SpaceTimeToy_Step (CCTK_PASS_FTOF)