diff options
Diffstat (limited to 'src/WaveToy.F')
| -rw-r--r-- | src/WaveToy.F | 192 |
1 files changed, 192 insertions, 0 deletions
diff --git a/src/WaveToy.F b/src/WaveToy.F new file mode 100644 index 0000000..f23ec94 --- /dev/null +++ b/src/WaveToy.F @@ -0,0 +1,192 @@ +#include "cctk.h" +#include "declare_parameters.h" +#include "declare_arguments.h" + + subroutine WaveToy_Boundaries(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_PARAMETERS + + INTEGER CCTK_Equals + + WaveToy_bound = "none" + + if(.NOT.CCTK_Equals(WaveToy_Bound, "periodic") .and. + & .NOT.CCTK_Equals(WaveToy_Bound,"none")) then + scalar_field(1,:,:) = 0.0 + scalar_field(:,1,:) = 0.0 + scalar_field(:,:,1) = 0.0 + + scalar_field(sh(1),:,:) = 0.0 + scalar_field(:,sh(2),:) = 0.0 + scalar_field(:,:,sh(3)) = 0.0 + end if + +c Should have a call to symmetry bounds here + + end subroutine WaveToy_boundaries + + + + subroutine WaveToy_evolution(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_PARAMETERS + +c Declare local variables + INTEGER :: i,j,k + INTEGER :: istart, jstart, kstart, iend, jend, kend + REAL :: dx,dy,dz + + dx = coarse_dx/levfac + dy = coarse_dy/levfac + dz = coarse_dz/levfac + + istart = 2 + jstart = 2 + kstart = 2 + + iend = sh(1)-1 + jend = sh(2)-1 + kend = sh(3)-1 + + do k = kstart, kend + do j = jstart, jend + do i = istart, iend + + scalar_field_tmp(i,j,k) = + 1 2.0*(1.0 - (dt**2)*(1.0/dx**2 + + 2 1.0/dy**2 +1.0/dz**2))*scalar_field(i,j,k) - + 3 scalar_field_old(i,j,k) + + 4 (dt**2) * + 5 ((scalar_field(i+1,j,k)+scalar_field(i-1,j,k))/dx**2 + 6 +(scalar_field(i,j+1,k)+scalar_field(i,j-1,k))/dy**2 + 7 +(scalar_field(i,j,k+1)+scalar_field(i,j,k-1))/dz**2) + + end do + end do + end do + + scalar_field_old = scalar_field + scalar_field = scalar_field_tmp + + call WaveToy_boundaries(CCTK_FARGUMENTS) + +c call CCTK_SyncGroup("scalarfields") + + end subroutine WaveToy_evolution + + + + subroutine WaveToy_init(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_PARAMETERS + + REAL :: omega, pi, kx, ky, kz, radius, amplitude + + integer CCTK_Equals + + pi = 4.0*atan(1.0) + +c Set parameters manually for now +c ------------------------------- + WaveToy_dtfac = 0.5 + WaveToy_kx = 2 + WaveToy_ky = 2 + WaveToy_kz = 2 + WaveToy_radius = 0.25 + WaveToy_amplitude = 1.0 + WaveToy_initialdata = "plane" + +c Calculate timestep + dt = WaveToy_dtfac*coarse_dx/levfac + time = 0 + +c Shorthand + kx = WaveToy_kx + ky = WaveToy_ky + kz = WaveToy_kz + radius = WaveToy_radius + amplitude = WaveToy_amplitude + + call WaveToy_boundaries(CCTK_FARGUMENTS) + +c call CCTK_SyncGroup("scalarfields") + + omega= sqrt(kx**2+ky**2+kz**2) + +c Initialise the scalar field. Can be plane waves, spherical waves +c or waves in a box. + + if (CCTK_Equals(WaveToy_InitialData,"plane")) then + scalar_field = amplitude* + & cos(kx*x+ky*y+kz*z+omega*time) + scalar_field_old = amplitude* + & cos(kx*x+ky*y+kz*z+omega*(time-dt)) + + else if (CCTK_Equals(WaveToy_initialdata,"spherical")) then + + scalar_field = amplitude*cos(kx*sqrt(x**2+y**2+z**2)+omega*time) + scalar_field_old = amplitude*cos(kx*sqrt(x**2+y**2+z**2)+omega*(time-dt)) + + else if (CCTK_Equals(WaveToy_InitialData, "box")) then + +c Only currently implemented for the "box" grid type. + +c if (CCTK_Equals(grid, "box")) then + +c Must have all the wavenumbers non-zero for this to work. + +c if (kx == 0) then +c write(*,*) "Box mode selected. Setting kx to 1" +c kx = 1 +c end if + +c if (ky == 0) then +c write(*,*) "Box mode selected. Setting ky to 1" +c ky = 1 +c end if + +c if (kz == 0) then +c write(*,*) "Box mode selected. Setting kz to 1" +c kz = 1 +c end if + +c Set up the initial data. +c Use kx,ky,kz as number of modes in each direction. + +c scalar_field = amplitude*sin(kx*(x-0.5)*pi)* +c $ sin(ky*(y-0.5)*pi)* +c $ sin(kz*(z-0.5)*pi)* +c $ cos(omega*time*pi) + +c scalar_field_old = amplitude*sin(kx*(x-0.5)*pi)* +c $ sin(ky*(y-0.5)*pi)* +c $ sin(kz*(z-0.5)*pi)* +c $ cos(omega*(time-dt)*pi) +c else +c write(*,*) "WaveToy box mode selected, but grid not set to box" +c STOP +c end if + + else if (CCTK_Equals(WaveToy_InitialData,"tsgaussian")) then + + scalar_field = exp(-amplitude*(sqrt(x**2+y**2+z**2)-radius)**2) + scalar_field_old = scalar_field + + else + + print *,"Unknown initialisation mode in WaveToy" +c CCTK_Warning(GH,0,"Unknown initialisation mode in WaveToy") +c CCTK_Stop(GH,0) + + end if + + end subroutine WaveToy_Init |