/*@@ @file WaveToy.F77 @date @author Tom Goodale @desc Evolution routines for the wave equation solver @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" /*@@ @routine WaveToyF77_Evolution @date @author Tom Goodale @desc Evolution for the wave equation @enddesc @calls CCTK_SyncGroup, wavetoy_boundaries @calledby @history @endhistory @@*/ subroutine WaveToyF77_Evolution(CCTK_ARGUMENTS) implicit none c Declare variables in argument list DECLARE_CCTK_ARGUMENTS INTEGER i,j,k,ierr INTEGER istart, jstart, kstart, iend, jend, kend CCTK_REAL dx,dy,dz,dt CCTK_REAL dx2,dy2,dz2,dt2 CCTK_REAL dx2i,dy2i,dz2i CCTK_REAL factor c Set up shorthands c ----------------- dx = CCTK_DELTA_SPACE(1) dy = CCTK_DELTA_SPACE(2) dz = CCTK_DELTA_SPACE(3) dt = CCTK_DELTA_TIME dx2 = dx*dx dy2 = dy*dy dz2 = dz*dz dt2 = dt*dt dx2i = 1.0/dx2 dy2i = 1.0/dy2 dz2i = 1.0/dz2 istart = 2 jstart = 2 kstart = 2 iend = cctk_lsh(1)-1 jend = cctk_lsh(2)-1 kend = cctk_lsh(3)-1 factor = 2*(1 - (dt2)*(dx2i + dy2i + dz2i)) c Do the evolution c ---------------- do k = kstart, kend do j = jstart, jend do i = istart, iend phi(i,j,k) = factor*phi_p(i,j,k) - & phi_p_p(i,j,k) + (dt2) * & ((phi_p(i+1,j,k)+phi_p(i-1,j,k))*dx2i & +(phi_p(i,j+1,k)+phi_p(i,j-1,k))*dy2i & +(phi_p(i,j,k+1)+phi_p(i,j,k-1))*dz2i) end do end do end do return end /*@@ @routine WaveToyF77_Boundaries @date @author Tom Goodale @desc Boundary conditions for the wave equation @enddesc @history @endhistory @@*/ subroutine WaveToyF77_Boundaries(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS c Local declarations CCTK_REAL zero,one CCTK_REAL finf integer npow integer ierr integer sw(3) ierr = -1 zero = 0.0 one = 1.0 npow = 1 finf = 1.0d0 c Set the stencil width c --------------------- sw(1)=1 sw(2)=1 sw(3)=1 c Apply the symmetry boundary conditions on any coordinate axes c ------------------------------------------------------------- call CartSymGN(ierr,cctkGH,"wavetoy::scalarevolve") c Apply the outer boundary conditions c ----------------------------------- if (CCTK_EQUALS(bound,"flat")) then call BndFlatVN(ierr,cctkGH,sw,"wavetoy::phi") else if (CCTK_EQUALS(bound,"zero")) then call BndScalarVN(ierr,cctkGH,zero,sw,"wavetoy::phi") else if (CCTK_Equals(bound,"radiation").eq.1) then call BndRadiativeVN(ierr,cctkGH,sw,zero,one,"wavetoy::phi", & "wavetoy::phi") else if (CCTK_Equals(bound,"robin").eq.1) then call BndRobinVN(ierr,cctkGH, sw, finf, npow,"wavetoy::phi") end if if (ierr < 0) then call CCTK_WARN(0,"Boundary conditions not applied - giving up!"); end if return end