diff options
| -rw-r--r-- | param.ccl | 1 | ||||
| -rw-r--r-- | src/GaugeWave.F | 124 |
2 files changed, 61 insertions, 64 deletions
@@ -440,6 +440,7 @@ KEYWORD GW_H "what function to use" { "sin" :: "1-a*sin(x)" "exp" :: "1-exp(x)" +"gau" :: "1-a*exp(-x**2)" }"sin" REAL GW_a "amplitude of the wave" diff --git a/src/GaugeWave.F b/src/GaugeWave.F index b11e6c6..78d6ddd 100644 --- a/src/GaugeWave.F +++ b/src/GaugeWave.F @@ -4,7 +4,7 @@ #define Pi (4 * atan(1.d0)) - + subroutine gauge_wave( $ x, y, z, t, $ gdtt, gdtx, gdty, gdtz, @@ -25,12 +25,17 @@ logical firstcall - CCTK_REAL a, o, H, half, d, fs + CCTK_REAL a, o, H, d, fs + CCTK_REAL zero,half,one data firstcall /.true./ - save firstcall, a, o, half, d, fs + save firstcall, a, o, half, d - +C Numbers. + + zero = 0.0d0 + half = 0.5d0 + one = 1.0d0 C Get parameters of the exact solution. @@ -38,25 +43,26 @@ C Get parameters of the exact solution. a = GW_a o = GW_omega - half = 1/2d0 -C d = GW_del / 8d0 / atan(1.d0) - d = GW_del * half / Pi + d = GW_del fs = GW_phase_shift firstcall = .false. end if -C How should the wave look like +C How should the wave look like. + if (CCTK_Equals(GW_H,"sin").eq.1) then - + + d = GW_del * half / Pi + if (GW_diagonal.ne.0) then - H = 1 - a * sin((x-y)/sqrt(2.d0)/d - o*t/d) + H = one - a * sin((x-y)/sqrt(2.d0)/d - o*t/d) else - H = 1 - a * sin((x-o*t)/d - fs ) + H = one - a * sin((x-o*t)/d - fs) end if @@ -64,65 +70,55 @@ C How should the wave look like H = exp(a*(x/d-o*t)) + elseif (CCTK_Equals(GW_H,"gau").eq.1) then + + H = one - a*dexp(-(x-t)**2/d**2) + end if -C write metric - - if (GW_diagonal.ne.0) then - gdxx = half + half *H - gdyy = half + half *H - gdxy = half - half *H - gdzz = 1 - gdzx = 0 - gdyz = 0 - - gdtt = - H - gdtx = 0 - gdty = 0 - gdtz = 0 - -C and upper metric - - guxx = half + half / H - guyy = half + half / H - guxy = half - half / H - guzz = 1 - guyz = 0 - guzx = 0 - - gutt = - 1/H - gutx = 0 - guty = 0 - gutz = 0 +C write metric. + + if (GW_diagonal .eq. 1) then + + gdxx = half * H + half + gdxy = -half* H + half + gdyy = half * H + half + + guxx = half / H + half + guxy =-half / H + half + guyy = half / H + half else - gdxx = H - gdyy = 1 - gdzz = 1 - gdxy = 0 - gdzx = 0 - gdyz = 0 - - gdtt = - H - gdtx = 0 - gdty = 0 - gdtz = 0 - -C upper metric - - guxx = 1/H - guyy = 1 - guzz = 1 - guxy = 0 - guyz = 0 - guzx = 0 - - gutt = - 1/H - gutx = 0 - guty = 0 - gutz = 0 + gdxx = H + gdxy = zero + gdyy = one + + guxx = one/H + guxy = zero + guyy = one + end if + + gdtt = - H + gdtx = zero + gdty = zero + gdtz = zero + + gdzx = zero + gdyz = zero + gdzz = one + +C and upper metric. + + gutt = - one/H + gutx = zero + guty = zero + gutz = zero + + guyz = zero + guzx = zero + guzz = one return end |