diff options
Diffstat (limited to 'CarpetExtra/IDScalarWave/src/InitialData.F77')
| -rw-r--r-- | CarpetExtra/IDScalarWave/src/InitialData.F77 | 121 |
1 files changed, 28 insertions, 93 deletions
diff --git a/CarpetExtra/IDScalarWave/src/InitialData.F77 b/CarpetExtra/IDScalarWave/src/InitialData.F77 index 0ee71e625..a6d74ae59 100644 --- a/CarpetExtra/IDScalarWave/src/InitialData.F77 +++ b/CarpetExtra/IDScalarWave/src/InitialData.F77 @@ -1,25 +1,23 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDScalarWave/src/InitialData.F77,v 1.14 2003/11/05 16:18:40 schnetter Exp $ /*@@ @file InitialData.F77 @date - @author Tom Goodale, Erik Schnetter + @author Tom Goodale @desc Initial data for the 3D Wave Equation @enddesc @@*/ #include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Functions.h" #include "cctk_Parameters.h" +#include "cctk_Arguments.h" /*@@ @routine IDScalarWave_InitialData @date - @author Tom Goodale, Erik Schnetter + @author Tom Goodale @desc Set up initial data for the wave equation @enddesc @@ -41,14 +39,11 @@ c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDScalarWave/src/Initi INTEGER i,j,k CCTK_REAL dt,omega, pi - CCTK_REAL t - CCTK_REAL x,y,z, r -c call CCTK_INFO ("IDScalarWave_InitialData") + print '("IDScalarWave_InitialData")' pi = 4.d0*atan(1.d0) - t = cctk_time dt = CCTK_DELTA_TIME omega = sqrt(kx**2+ky**2+kz**2) @@ -59,18 +54,12 @@ c call CCTK_INFO ("IDScalarWave_InitialData") 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) - phi(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) phi_p(i,j,k) = amplitude - $ * cos((kx*x + ky*y + kz*z + omega*(cctk_time - dt)) * pi) - - phi_p_p(i,j,k) = amplitude - $ * cos((kx*x + ky*y + kz*z + omega*(cctk_time - 2*dt)) * pi) + $ * cos(kx*cart3d_x(i,j,k) + ky*cart3d_y(i,j,k) + $ + kz*cart3d_z(i,j,k) + omega*(cctk_time - dt)) end do end do @@ -81,41 +70,14 @@ c call CCTK_INFO ("IDScalarWave_InitialData") do k=1, cctk_lsh(3) do j=1, cctk_lsh(2) do i=1, cctk_lsh(1) - - r = spher3d_r(i,j,k) - - if (r .ne. 0) then - phi(i,j,k) = amplitude/2 * (r - radius - t) / r - $ * exp(- (r - radius - t)**2 / sigma**2) - $ + amplitude/2 * (r + radius + t) / r - $ * exp(- (r + radius + t)**2 / sigma**2) - - phi_p(i,j,k) = amplitude/2 * (r - radius - t + dt) / r - $ * exp(- (r - radius - t + dt)**2 / sigma**2) - $ + amplitude/2 * (r + radius + t - dt) / r - $ * exp(- (r + radius + t - dt)**2 / sigma**2) - - phi_p_p(i,j,k) = amplitude/2 * (r - radius - t + 2*dt) / r - $ * exp(- (r - radius - t + 2*dt)**2 / sigma**2) - $ + amplitude/2 * (r + radius + t - 2*dt) / r - $ * exp(- (r + radius + t - 2*dt)**2 / sigma**2) - - else - - phi(i,j,k) = amplitude - $ * (1 - 2 * (radius + t)**2 / sigma**2) - $ * exp(- (radius + t)**2 / sigma**2) - - phi_p(i,j,k) = amplitude - $ * (1 - 2 * (radius + t - dt)**2 / sigma**2) - $ * exp(- (radius + t - dt)**2 / sigma**2) - - phi_p_p(i,j,k) = amplitude - $ * (1 - 2 * (radius + t - 2*dt)**2 / sigma**2) - $ * exp(- (radius + t - 2*dt)**2 / sigma**2) - - end if - + + phi(i,j,k) = amplitude + $ * exp(- (sqrt(cart3d_x(i,j,k)**2 + cart3d_y(i,j,k)**2 + $ + cart3d_z(i,j,k)**2) - radius)**2 / sigma**2) + phi_p(i,j,k) = amplitude + $ * exp(- (sqrt(cart3d_x(i,j,k)**2 + cart3d_y(i,j,k)**2 + $ + cart3d_z(i,j,k)**2) - radius - dt)**2 / sigma**2) + end do end do end do @@ -128,59 +90,30 @@ 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) - phi(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) - phi_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) - - phi_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) + 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) 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) - - x = cart3d_x(i,j,k) - y = cart3d_y(i,j,k) - z = cart3d_z(i,j,k) - - phi(i,j,k) = 1 / sqrt((x-cx)**2 + (y-cy)**2 + (z-cz)**2) - 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 + phi(i,j,k) = 0.0d0 + phi_p(i,j,k) = 0.0d0 end do end do @@ -189,3 +122,5 @@ c Use kx,ky,kz as number of modes in each direction. end if end + + |