diff options
Diffstat (limited to 'src/InitialData.F')
-rw-r--r-- | src/InitialData.F | 96 |
1 files changed, 96 insertions, 0 deletions
diff --git a/src/InitialData.F b/src/InitialData.F new file mode 100644 index 0000000..1f12a5f --- /dev/null +++ b/src/InitialData.F @@ -0,0 +1,96 @@ + /*@@ + @file InitialData.F + @date + @author Tom Goodale + @desc + Initial data for the 3D Wave Equation + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + + + /*@@ + @routine WaveToyF90_InitialData + @date + @author Tom Goodale + @desc + Set up initial data for the wave equation + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + + subroutine WaveToyF90_InitialData(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + INTEGER CCTK_Equals + + INTEGER :: i + CCTK_REAL :: dt,omega, pi + CCTK_REAL :: min_delta,dx,dy,dz + + pi = 4.0*atan(1.0) + +c Grid spacing shortcuts +c ---------------------- + dx = cctk_delta_space(1) + dy = cctk_delta_space(2) + dz = cctk_delta_space(3) + +c Calculate timestep +c ------------------ + min_delta = min(dx,dy,dz) + cctk_delta_time = dtfac*min_delta + dt = cctk_delta_time + omega = sqrt(kx**2+ky**2+kz**2) + + if (CCTK_Equals(initial_data,"plane")==1) then + + phi = amplitude*cos(kx*x+ky*y+kz*z+omega*cctk_time) + phi_old = amplitude*cos(kx*x+ky*y+kz*z+omega*(cctk_time-dt)) + + else if (CCTK_Equals(initial_data,"gaussian")==1) then + + call CCTK_INFO("Gaussian initial data for Wave Equation"); + phi = amplitude*exp( -(sqrt(x**2+y**2+z**2)-radius)**2/sigma**2) + phi_old = amplitude*exp( -(sqrt(x**2+y**2+z**2)-radius-dt)**2/sigma**2) + + else if (CCTK_Equals(initial_data, "box")==1) then + +c Use kx,ky,kz as number of modes in each direction. + + phi = amplitude*sin(kx*(x-0.5)*pi)* + $ sin(ky*(y-0.5)*pi)* + $ sin(kz*(z-0.5)*pi)* + $ cos(omega*cctk_time*pi) + + phi_old = amplitude*sin(kx*(x-0.5)*pi)* + $ sin(ky*(y-0.5)*pi)* + $ sin(kz*(z-0.5)*pi)* + $ cos(omega*(cctk_time-dt)*pi) + + end if + +c Apply symmetry boundary conditions +c ---------------------------------- + call ApplySymmetry(cctkGH,"wavetoy::scalarevolve") + +c Synchronise +c ----------- + call CCTK_SyncGroup(cctkGH,"wavetoy::scalarevolve") + + end subroutine WaveToyF90_InitialData + + |