Added Gowdy periodic wave solution.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/Exact/trunk@145 e296648e-0e4f-0410-bd07-d597d9acff87

7 files changed, 262 insertions, 1 deletions

diff --git a/param.ccl b/param.ccl
index 3e79e97..2f26830 100644
--- a/param.ccl
+++ b/param.ccl
@@ -188,6 +188,7 @@ KEYWORD exact_model "The exact solution/coordinates used in thorn exact"
 "Kasner-like"          :: "Kasner-like spacetime"
 "Kasner-axisymmetric"  :: "axisymmetric Kasner spacetime"
 "Kasner-generalized"   :: "generalized Kasner spacetime"
+"Gowdy-wave"	       :: "Gowdy spacetime with polarized wave on a torus"
 "Milne"                :: "Milne spacetime for pre-big-bang cosmology"
 #
 # miscelaneous spacetimes
@@ -707,6 +708,17 @@ REAL Kasner_generalized__p2 "Kasner-generalized: y exponent parameter"
 ################################################################################
 
 #
+# parameters for Gowdy polarized wave spacetime
+#
+
+REAL Gowdy_wave__amplitude "Gowdy-wave: amplitude parameter"
+{
+*:* :: "any real number"
+} 0.0
+
+################################################################################
+
+#
 # parameters for Milne spacetime
 #
 
diff --git a/src/decode_pars.F77 b/src/decode_pars.F77
index e81fd8a..1368f11 100644
--- a/src/decode_pars.F77
+++ b/src/decode_pars.F77
@@ -101,6 +101,8 @@ c cosmological spacetimes
 	     decoded_exact_model = EXACT__Kasner_axisymmetric
 	elseif (CCTK_Equals(exact_model, "Kasner-generalized") .ne. 0) then
 	     decoded_exact_model = EXACT__Kasner_generalized
+	elseif (CCTK_Equals(exact_model, "Gowdy-wave") .ne. 0) then
+	     decoded_exact_model = EXACT__Gowdy_wave
 	elseif (CCTK_Equals(exact_model, "Milne") .ne. 0) then
 	     decoded_exact_model = EXACT__Milne
 
diff --git a/src/include/Scalar_CalcTmunu.inc b/src/include/Scalar_CalcTmunu.inc
index 5bb7cbf..4f5f240 100644
--- a/src/include/Scalar_CalcTmunu.inc
+++ b/src/include/Scalar_CalcTmunu.inc
@@ -287,6 +287,11 @@ c
         Tyz  = Tyz
 
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+       elseif (decoded_exact_model .eq. EXACT__Gowdy_wave) then
+c        no stress-energy tensor in this model
+
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
        elseif (decoded_exact_model .eq. EXACT__Milne) then
 c        no stress-energy tensor in this model
diff --git a/src/include/param_defs.inc b/src/include/param_defs.inc
index 2b4bbd2..df25a06 100644
--- a/src/include/param_defs.inc
+++ b/src/include/param_defs.inc
@@ -58,7 +58,8 @@ c cosmological spacetimes
 #define EXACT__Kasner_like		58
 #define EXACT__Kasner_axisymmetric	59
 #define EXACT__Kasner_generalized	60
-#define EXACT__Milne			61
+#define EXACT__Gowdy_wave		61
+#define EXACT__Milne			62
 
 c miscelaneous spacetimes
 #define EXACT__boost_rotation_symmetric	80
diff --git a/src/metric.F77 b/src/metric.F77
index 23797f7..c547b8d 100644
--- a/src/metric.F77
+++ b/src/metric.F77
@@ -243,6 +243,14 @@ c
      $        gutt, gutx, guty, gutz,
      $        guxx, guyy, guzz, guxy, guyz, guxz)
 
+      elseif (decoded_exact_model .eq. EXACT__Gowdy_wave) then
+         call Exact__Gowdy_wave(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
       elseif (decoded_exact_model .eq. EXACT__Milne) then
          call Exact__Milne(
      $        x, y, z, t,
diff --git a/src/metrics/Gowdy.F77 b/src/metrics/Gowdy.F77
new file mode 100644
index 0000000..06370bf
--- /dev/null
+++ b/src/metrics/Gowdy.F77
@@ -0,0 +1,232 @@
+C @@
+C  @file      Gowdy.F77
+C  @date      December 2002
+C  @author    Denis Pollney
+C  @desc
+C             Cosmological Gowdy metric for a polarized wave in an
+C	      expanding universe. See
+C	        "Stable 3-level leapfrog integration in numerical
+C		relativity" New, K, Watt, K, Misner, C and Centrella,
+C		J, PRD 58, 064022.
+C  @desc
+C  @version   $Header$
+C @@
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+      subroutine Exact__Gowdy_wave(
+     $     x, y, z, t,
+     $     gdtt, gdtx, gdty, gdtz, 
+     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
+     $     gutt, gutx, guty, gutz, 
+     $     guxx, guyy, guzz, guxy, guyz, guzx)
+
+      implicit none
+
+      DECLARE_CCTK_PARAMETERS
+
+      CCTK_REAL x, y, z, t
+      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
+     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
+     $       gutt, gutx, guty, gutz, 
+     $       guxx, guyy, guzz, guxy, guyz, guzx
+
+      logical firstcall
+    
+      CCTK_REAL amp
+      CCTK_REAL Bessel_J0, Bessel_J1, Bessel_J0_t, Bessel_J1_t
+      CCTK_REAL PI, twoPI, cosz, eP, lambda
+      CCTK_REAL d1, d2, d3, d4, d5, d6
+
+      data firstcall /.true./
+      save firstcall, amp, PI, twoPI, Bessel_J0, Bessel_J1
+
+      if (firstcall) then
+          amp = Gowdy_wave__amplitude
+	  PI = acos(-1.d0)
+	  twoPI = 2.d0*PI
+	  call jy01a(twoPI, Bessel_J0, d1, Bessel_J1, d2, d3, d4, d5, d6)
+          firstcall = .false.
+      end if
+
+      if (t.eq.0.d0) then
+	call CCTK_WARN(0, "The Gowdy metric is singular for t=0. You may need to set cactus::cctk_initial_time > 0.")
+      end if
+
+      call jy01a(twoPI*t, Bessel_J0_t, d1, Bessel_J1_t, d2, d3, d4, d5, d6)
+      
+      cosz = cos(twoPI*z)
+      eP = exp(Bessel_J0_t * cosz)
+
+      lambda = amp * (-twoPI * t * Bessel_J0_t * Bessel_J1_t * cosz**2
+     +  + twoPI * PI * t**2 * (Bessel_J0_t**2 + Bessel_J1_t**2)
+     +  - 0.5d0 * (twoPI**2 * (Bessel_J0**2 + Bessel_J1**2)
+     +             - twoPI * Bessel_J0 * Bessel_J1))
+c
+c     lower metric components.
+c
+      gdtt = -exp(0.5d0*lambda)/sqrt(t)
+      gdtx = 0.d0 
+      gdty = 0.d0
+      gdtz = 0.d0
+      gdxx = t * eP
+      gdyy = t / eP
+      gdzz = -gdtt
+      gdxy = 0.d0
+      gdyz = 0.d0
+      gdzx = 0.d0
+
+c
+c     upper metric components.
+c
+      gutt = 1.d0 / gdtt
+      gutx = 0.d0
+      guty = 0.d0
+      gutz = 0.d0
+      guxx = 1.d0 / gdxx
+      guyy = 1.d0 / gdyy
+      guzz = 1.d0 / gdzz
+      guxy = 0.d0
+      guyz = 0.d0
+      guzx = 0.d0
+
+      return
+      end
+
+C @@
+C  @routine   Bessel.F77
+C  @date      December 2002
+C  @author    Denis Pollney
+C  @desc
+C             Compute bessel functions of 0th and 1st order.
+C             This routine was downloaded from:
+C               http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
+C
+C       =======================================================
+C       Purpose: Compute Bessel functions J0(x), J1(x), Y0(x),
+C                Y1(x), and their derivatives
+C       Input :  x   --- Argument of Jn(x) & Yn(x) ( x ? 0 )
+C       Output:  BJ0 --- J0(x)
+C                DJ0 --- dJ0(x)
+C                BJ1 --- J1(x)
+C                DJ1 --- dJ1(x)
+C                BY0 --- Y0(x)
+C                DY0 --- dY0(x)
+C                BY1 --- Y1(x)
+C                DY1 --- dY1(x)
+C       =======================================================
+C
+C  @desc
+C  @version   $Header$
+C @@
+
+        SUBROUTINE JY01A(X,BJ0,DJ0,BJ1,DJ1,BY0,DY0,BY1,DY1)
+        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+        DIMENSION A(12),B(12),A1(12),B1(12)
+        PI=3.141592653589793D0
+        RP2=0.63661977236758D0
+        X2=X*X
+        IF (X.EQ.0.0D0) THEN
+           BJ0=1.0D0
+           BJ1=0.0D0
+           DJ0=0.0D0
+           DJ1=0.5D0
+           BY0=-1.0D+300
+           BY1=-1.0D+300
+           DY0=1.0D+300
+           DY1=1.0D+300
+           RETURN
+        ENDIF
+        IF (X.LE.12.0D0) THEN
+           BJ0=1.0D0
+           R=1.0D0
+           DO 5 K=1,30
+              R=-0.25D0*R*X2/(K*K)
+              BJ0=BJ0+R
+              IF (DABS(R).LT.DABS(BJ0)*1.0D-15) GO TO 10
+5          CONTINUE
+10         BJ1=1.0D0
+           R=1.0D0
+           DO 15 K=1,30
+              R=-0.25D0*R*X2/(K*(K+1.0D0))
+              BJ1=BJ1+R
+              IF (DABS(R).LT.DABS(BJ1)*1.0D-15) GO TO 20
+15         CONTINUE
+20         BJ1=0.5D0*X*BJ1
+           EC=DLOG(X/2.0D0)+0.5772156649015329D0
+           CS0=0.0D0
+           W0=0.0D0
+           R0=1.0D0
+           DO 25 K=1,30
+              W0=W0+1.0D0/K
+              R0=-0.25D0*R0/(K*K)*X2
+              R=R0*W0
+              CS0=CS0+R
+              IF (DABS(R).LT.DABS(CS0)*1.0D-15) GO TO 30
+25         CONTINUE
+30         BY0=RP2*(EC*BJ0-CS0)
+           CS1=1.0D0
+           W1=0.0D0
+           R1=1.0D0
+           DO 35 K=1,30
+              W1=W1+1.0D0/K
+              R1=-0.25D0*R1/(K*(K+1))*X2
+              R=R1*(2.0D0*W1+1.0D0/(K+1.0D0))
+              CS1=CS1+R
+              IF (DABS(R).LT.DABS(CS1)*1.0D-15) GO TO 40
+35         CONTINUE
+40         BY1=RP2*(EC*BJ1-1.0D0/X-0.25D0*X*CS1)
+        ELSE
+           DATA A/-.7031250000000000D-01,.1121520996093750D+00,
+     &            -.5725014209747314D+00,.6074042001273483D+01,
+     &            -.1100171402692467D+03,.3038090510922384D+04,
+     &            -.1188384262567832D+06,.6252951493434797D+07,
+     &            -.4259392165047669D+09,.3646840080706556D+11,
+     &            -.3833534661393944D+13,.4854014686852901D+15/
+           DATA B/ .7324218750000000D-01,-.2271080017089844D+00,
+     &             .1727727502584457D+01,-.2438052969955606D+02,
+     &             .5513358961220206D+03,-.1825775547429318D+05,
+     &             .8328593040162893D+06,-.5006958953198893D+08,
+     &             .3836255180230433D+10,-.3649010818849833D+12,
+     &             .4218971570284096D+14,-.5827244631566907D+16/
+           DATA A1/.1171875000000000D+00,-.1441955566406250D+00,
+     &             .6765925884246826D+00,-.6883914268109947D+01,
+     &             .1215978918765359D+03,-.3302272294480852D+04,
+     &             .1276412726461746D+06,-.6656367718817688D+07,
+     &             .4502786003050393D+09,-.3833857520742790D+11,
+     &             .4011838599133198D+13,-.5060568503314727D+15/
+           DATA B1/-.1025390625000000D+00,.2775764465332031D+00,
+     &             -.1993531733751297D+01,.2724882731126854D+02,
+     &             -.6038440767050702D+03,.1971837591223663D+05,
+     &             -.8902978767070678D+06,.5310411010968522D+08,
+     &             -.4043620325107754D+10,.3827011346598605D+12,
+     &             -.4406481417852278D+14,.6065091351222699D+16/
+           K0=12
+           IF (X.GE.35.0) K0=10
+           IF (X.GE.50.0) K0=8
+           T1=X-0.25D0*PI
+           P0=1.0D0
+           Q0=-0.125D0/X
+           DO 45 K=1,K0
+              P0=P0+A(K)*X**(-2*K)
+45            Q0=Q0+B(K)*X**(-2*K-1)
+           CU=DSQRT(RP2/X)
+           BJ0=CU*(P0*DCOS(T1)-Q0*DSIN(T1))
+           BY0=CU*(P0*DSIN(T1)+Q0*DCOS(T1))
+           T2=X-0.75D0*PI
+           P1=1.0D0
+           Q1=0.375D0/X
+           DO 50 K=1,K0
+              P1=P1+A1(K)*X**(-2*K)
+50            Q1=Q1+B1(K)*X**(-2*K-1)
+           CU=DSQRT(RP2/X)
+           BJ1=CU*(P1*DCOS(T2)-Q1*DSIN(T2))
+           BY1=CU*(P1*DSIN(T2)+Q1*DCOS(T2))
+        ENDIF
+        DJ0=-BJ1
+        DJ1=BJ0-BJ1/X
+        DY0=-BY1
+        DY1=BY0-BY1/X
+        RETURN
+        END
diff --git a/src/metrics/make.code.defn b/src/metrics/make.code.defn
index adabc77..9c33332 100644
--- a/src/metrics/make.code.defn
+++ b/src/metrics/make.code.defn
@@ -31,6 +31,7 @@ SRCS = Minkowski.F77			\
        Kasner_like.F77			\
        Kasner_axisymmetric.F77		\
        Kasner_generalized.F77		\
+       Gowdy.F77			\
        Milne.F77			\
 					\
        boost_rotation_symmetric.F77	\