Initial revision

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

2 files changed, 288 insertions, 0 deletions

diff --git a/CarpetExtra/IDFOScalarWave/src/InitialData.F77 b/CarpetExtra/IDFOScalarWave/src/InitialData.F77
new file mode 100644
index 000000000..aa361f8a1
--- /dev/null
+++ b/CarpetExtra/IDFOScalarWave/src/InitialData.F77
@@ -0,0 +1,279 @@
+c     -*-Fortran-*-
+c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDFOScalarWave/src/InitialData.F77,v 1.8 2003/11/05 16:18:40 schnetter Exp $
+
+ /*@@
+   @file      InitialData.F77
+   @date      
+   @author    Scott Hawley, using code from Tom Goodale, Erik Schnetter
+   @desc 
+              Initial data for the 3D Wave Equation
+   @enddesc 
+ @@*/
+
+#include "cctk.h" 
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+
+ /*@@
+   @routine    IDFOScalarWave_InitialData
+   @date       
+   @author     Scott Hawley, using code from Tom Goodale, Erik Schnetter
+   @desc 
+               Set up initial data for the wave equation
+   @enddesc 
+   @calls      
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+      subroutine IDFOScalarWave_InitialData (CCTK_ARGUMENTS)
+
+      implicit none
+
+      DECLARE_CCTK_ARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
+      INTEGER    i,j,k
+      CCTK_REAL  dt,omega, cpi
+      CCTK_REAL  ri3
+      
+c     call CCTK_INFO ("IDFOScalarWave_InitialData")
+      
+      cpi = 4.d0*atan(1.d0)
+
+      dt = CCTK_DELTA_TIME
+
+      omega = sqrt(kx**2+ky**2+kz**2)
+
+      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)
+
+                  phi(i,j,k) = amplitude
+     $                 * cos((kx*x(i,j,k) + ky*y(i,j,k) + kz*z(i,j,k) 
+     $                 + omega*cctk_time) * cpi)
+
+                  phi_p(i,j,k) = amplitude
+     $              * cos((kx*x(i,j,k) + ky*y(i,j,k) + kz*z(i,j,k) 
+     $              + omega*(cctk_time - dt)) * cpi)
+
+                  phi_p_p(i,j,k) = amplitude
+     $             * cos((kx*x(i,j,k) + ky*y(i,j,k) + kz*z(i,j,k) 
+     $              + omega*(cctk_time - 2*dt)) * cpi)
+                  
+               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)
+
+                  phi(i,j,k) = amplitude 
+     $                 * exp(- (x(i,j,k) - radius)**2 / sigma**2)
+     $                 * exp(- (y(i,j,k) - radius)**2 / sigma**2)
+     $                 * exp(- (z(i,j,k) - radius)**2 / sigma**2)
+
+                  pi(i,j,k) = 0.0
+                  phix(i,j,k) = phi(i,j,k)*  (-2) * (x(i,j,k) - radius) 
+     $              / sigma**2
+                  phiy(i,j,k) = phi(i,j,k)*  (-2) * (y(i,j,k) - radius) 
+     $              / sigma**2
+                  phiz(i,j,k) = phi(i,j,k)*  (-2) * (z(i,j,k) - radius) 
+     $              / sigma**2
+
+                  pi_p(i,j,k) = amplitude *
+     &                    exp( -(x(i,j,k) - radius - dt)**2 / sigma**2 )* 
+     &                    exp( -(y(i,j,k) - radius - dt)**2 / sigma**2 )* 
+     &                    exp( -(z(i,j,k) - radius - dt)**2 / sigma**2 )* 
+     &               ((x(i,j,k)-radius-dt) + (y(i,j,k)-radius-dt) 
+     $              + (z(i,j,k)-radius-dt))/sigma**2
+     &                        - amplitude * 
+     &                    exp( -(x(i,j,k) - radius + dt)**2 / sigma**2 )* 
+     &                    exp( -(y(i,j,k) - radius + dt)**2 / sigma**2 )* 
+     &                    exp( -(z(i,j,k) - radius + dt)**2 / sigma**2 )* 
+     &               ((x(i,j,k)-radius+dt) + (y(i,j,k)-radius+dt) 
+     &              + (z(i,j,k)-radius+dt))/sigma**2
+
+                  pi_p_p(i,j,k) = amplitude * 
+     &                    exp( -(x(i,j,k) - radius - 2*dt)**2 / sigma**2 )* 
+     &                    exp( -(y(i,j,k) - radius - 2*dt)**2 / sigma**2 )* 
+     &                    exp( -(z(i,j,k) - radius - 2*dt)**2 / sigma**2 )* 
+     &          ((x(i,j,k)-radius-2*dt) + (y(i,j,k)-radius-2*dt) 
+     &              + (z(i,j,k)-radius-2*dt))/sigma**2
+     &                        - amplitude * 
+     &                    exp( -(x(i,j,k) - radius + 2*dt)**2 / sigma**2 )* 
+     &                    exp( -(y(i,j,k) - radius + 2*dt)**2 / sigma**2 )* 
+     &                    exp( -(z(i,j,k) - radius + 2*dt)**2 / sigma**2 )* 
+     &          ((x(i,j,k)-radius+2*dt) + (y(i,j,k)-radius+2*dt) 
+     &              + (z(i,j,k)-radius+2*dt))/sigma**2
+
+                  phix_p(i,j,k) = - amplitude * (x(i,j,k) - radius - dt) 
+     &              / sigma**2 
+     &                 *  exp( -(x(i,j,k) - radius - dt)**2 / sigma**2 )* 
+     &                    exp( -(y(i,j,k) - radius - dt)**2 / sigma**2 )* 
+     &                    exp( -(z(i,j,k) - radius - dt)**2 / sigma**2 )
+     &                            - amplitude * (x(i,j,k) - radius + dt) 
+     &              / sigma**2
+     &                 *  exp( -(x(i,j,k) - radius + dt)**2 / sigma**2 )* 
+     &                    exp( -(y(i,j,k) - radius + dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius + dt)**2 / sigma**2 )
+     
+                  phix_p_p(i,j,k) = - amplitude * (x(i,j,k) - radius - 2*dt) 
+     &              / sigma**2
+     &                 *  exp( -(x(i,j,k) - radius - 2*dt)**2 / sigma**2 )* 
+     &                    exp( -(y(i,j,k) - radius - 2*dt)**2 / sigma**2 )* 
+     &                    exp( -(z(i,j,k) - radius - 2*dt)**2 / sigma**2 )
+     &                            - amplitude * (x(i,j,k) - radius + 2*dt) 
+     &              / sigma**2
+     &                 *  exp( -(x(i,j,k) - radius + 2*dt)**2 / sigma**2 )* 
+     &                    exp( -(y(i,j,k) - radius + 2*dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius + 2*dt)**2 / sigma**2 )
+
+                  phiy_p(i,j,k) = - amplitude * (y(i,j,k) - radius - dt) 
+     &              / sigma**2
+     &                 *  exp( -(x(i,j,k) - radius - dt)**2 / sigma**2 )*
+     &                    exp( -(y(i,j,k) - radius - dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius - dt)**2 / sigma**2 ) 
+     &                            - amplitude * (y(i,j,k) - radius + dt) 
+     &              / sigma**2
+     &                 *  exp( -(x(i,j,k) - radius + dt)**2 / sigma**2 )*
+     &                    exp( -(y(i,j,k) - radius + dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius + dt)**2 / sigma**2 ) 
+    
+                  phiy_p_p(i,j,k) = - amplitude * (y(i,j,k) - radius - 2*dt) 
+     &              / sigma**2
+     &                 *  exp( -(x(i,j,k) - radius - 2*dt)**2 / sigma**2 )*
+     &                    exp( -(y(i,j,k) - radius - 2*dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius - 2*dt)**2 / sigma**2 )
+     &                            - amplitude * (y(i,j,k) - radius + 2*dt) 
+     &              / sigma**2
+     &                 *  exp( -(x(i,j,k) - radius + 2*dt)**2 / sigma**2 )*
+     &                    exp( -(y(i,j,k) - radius + 2*dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius + 2*dt)**2 / sigma**2 )
+
+                  phiz_p(i,j,k) = - amplitude * (z(i,j,k) - radius - dt) 
+     &              / sigma**2
+     &                 *  exp( -(x(i,j,k) - radius - dt)**2 / sigma**2 )*
+     &                    exp( -(y(i,j,k) - radius - dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius - dt)**2 / sigma**2 ) 
+     &                            - amplitude * (z(i,j,k) - radius + dt) 
+     &              / sigma**2
+     &                 *  exp( -(x(i,j,k) - radius + dt)**2 / sigma**2 )*
+     &                    exp( -(y(i,j,k) - radius + dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius + dt)**2 / sigma**2 ) 
+    
+                  phiz_p_p(i,j,k) = - amplitude * (z(i,j,k) - radius - 2*dt) 
+     &              / sigma**2
+     &                 *  exp( -(x(i,j,k) - radius - 2*dt)**2 / sigma**2 )*
+     &                    exp( -(y(i,j,k) - radius - 2*dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius - 2*dt)**2 / sigma**2 )
+     &                            - amplitude * (z(i,j,k) - radius + 2*dt) 
+     &              / sigma**2
+     &                 *  exp( -(x(i,j,k) - radius + 2*dt)**2 / sigma**2 )*
+     &                    exp( -(y(i,j,k) - radius + 2*dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius + 2*dt)**2 / sigma**2 )
+
+                  phi_p(i,j,k) = amplitude / 2 *
+     &                    exp( -(x(i,j,k) - radius - dt)**2 / sigma**2 )*
+     &                    exp( -(y(i,j,k) - radius - dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius - dt)**2 / sigma**2 )
+     &                          + amplitude / 2 *
+     &                    exp( -(x(i,j,k) - radius + dt)**2 / sigma**2 )*
+     &                    exp( -(y(i,j,k) - radius + dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius + dt)**2 / sigma**2 )
+
+                  phi_p_p(i,j,k) = amplitude / 2 *
+     &                    exp( -(x(i,j,k) - radius - 2*dt)**2 / sigma**2 )*
+     &                    exp( -(y(i,j,k) - radius - 2*dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius - 2*dt)**2 / sigma**2 )
+     &                          + amplitude / 2 *
+     &                    exp( -(x(i,j,k) - radius + 2*dt)**2 / sigma**2 )*
+     &                    exp( -(y(i,j,k) - radius + 2*dt)**2 / sigma**2 )*
+     &                    exp( -(z(i,j,k) - radius + 2*dt)**2 / sigma**2 )
+
+               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)
+
+                  phi(i,j,k) = amplitude
+     $                 * sin(kx * (x(i,j,k) - 0.5d0) * cpi)
+     $                 * sin(ky * (y(i,j,k) - 0.5d0) * cpi)
+     $                 * sin(kz * (z(i,j,k) - 0.5d0) * cpi)
+                  
+                  phi_p(i,j,k) = phi(i,j,k)
+     $                 * cos(omega * (cctk_time - dt) * cpi)
+                  
+                  phi_p_p(i,j,k) = phi(i,j,k)
+     $                 * cos(omega * (cctk_time - 2*dt) * cpi)
+
+                  phi(i,j,k) = phi(i,j,k) * cos(omega * cctk_time * cpi)
+
+               end do
+            end do
+         end do
+
+      else if (CCTK_EQUALS(initial_data, "1/r")) then
+  
+         do k=1,cctk_lsh(3)
+            do j=1,cctk_lsh(2)
+               do i=1,cctk_lsh(1)
+                  
+                  pi(i,j,k)  = 0.0
+                  phi(i,j,k) = 1 / sqrt(x(i,j,k)**2 + y(i,j,k)**2 
+     &                     + z(i,j,k)**2)
+                  ri3 = phi(i,j,k)**3
+                  phix(i,j,k) = - x(i,j,k) * ri3
+                  phiy(i,j,k) = - y(i,j,k) * ri3
+                  phiz(i,j,k) = - z(i,j,k) * ri3
+
+                  pi_p(i,j,k)     = pi(i,j,k)
+                  pi_p_p(i,j,k)   = pi(i,j,k)
+                  phix_p(i,j,k)   = phix(i,j,k)
+                  phix_p_p(i,j,k) = phix(i,j,k)
+                  phiy_p(i,j,k)   = phiy(i,j,k)
+                  phiy_p_p(i,j,k) = phiy(i,j,k)
+                  phiz_p(i,j,k)   = phiz(i,j,k)
+                  phiz_p_p(i,j,k) = phiz(i,j,k)
+                  phi_p(i,j,k)    = phi(i,j,k)
+                  phi_p_p(i,j,k)  = phi(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)
+                  
+                  phi(i,j,k)     = 0.0d0
+                  phi_p(i,j,k)   = 0.0d0
+                  phi_p_p(i,j,k) = 0.0d0
+
+               end do
+            end do
+         end do 
+         
+      end if
+      
+      end 
diff --git a/CarpetExtra/IDFOScalarWave/src/make.code.defn b/CarpetExtra/IDFOScalarWave/src/make.code.defn
new file mode 100644
index 000000000..3d9912022
--- /dev/null
+++ b/CarpetExtra/IDFOScalarWave/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn IDScalarWave   -*-Makefile-*-
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDFOScalarWave/src/make.code.defn,v 1.2 2003/06/27 16:26:11 schnetter Exp $
+
+# Source files in this directory
+SRCS = InitialData.F77
+
+# Subdirectories containing source files
+SUBDIRS = 
+