diff options
Diffstat (limited to 'src/WaveToy.c')
-rw-r--r-- | src/WaveToy.c | 157 |
1 files changed, 157 insertions, 0 deletions
diff --git a/src/WaveToy.c b/src/WaveToy.c new file mode 100644 index 0000000..5635056 --- /dev/null +++ b/src/WaveToy.c @@ -0,0 +1,157 @@ + /*@@ + @file WaveToy.c + @date + @author Tom Goodale + @desc + Evolution routines for the wave equation solver + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + +#include "cctk_Flesh.h" +#include "cctk_Misc.h" +#include "cctk_Comm.h" + +#include "CactusBase/Boundary/src/Boundary.h" + + /*@@ + @routine WaveToyC_Boundaries + @date + @author Tom Goodale + @desc + Boundary conditions for the wave equation + @enddesc + @calls ApplyFlatBC,ApplyRadiativeBC + @calledby + @history + + @endhistory + +@@*/ + +void WaveToyC_Boundaries(CCTK_CARGUMENTS) +{ + DECLARE_CCTK_CARGUMENTS + DECLARE_CCTK_PARAMETERS + + int sw[3]; + + /* Set the stencil width */ + sw[0]=1; + sw[1]=1; + sw[2]=1; + + ApplySymmetry(cctkGH,"wavetoy::scalarevolve"); + + if (CCTK_EQUALS(bound,"flat")) + { + ApplyFlatBC(cctkGH,sw,"wavetoy::phi"); + } + else if (CCTK_Equals(bound,"radiation")) + { + ApplyRadiativeBC(cctkGH,1,sw,"wavetoy::phi","wavetoy::phi_old"); + } + +} + + + /*@@ + @routine WaveToyC_Evolution + @date + @author Tom Goodale + @desc + Evolution for the wave equation + @enddesc + @calls CCTK_SyncGroup, wavetoy_boundaries + @calledby + @history + + @endhistory + +@@*/ + +void WaveToyC_Evolution(CCTK_CARGUMENTS) +{ + DECLARE_CCTK_CARGUMENTS + DECLARE_CCTK_PARAMETERS + + int i,j,k; + int istart, jstart, kstart, iend, jend, kend; + CCTK_REAL dx,dy,dz,dt,dx2,dy2,dz2,dt2; + + /* Set up shorthands */ + dx = cctk_delta_space[0]; + dy = cctk_delta_space[1]; + dz = cctk_delta_space[2]; + dt = cctk_delta_time; + + dx2=dx*dx; + dy2=dy*dy; + dz2=dz*dz; + dt2=dt*dt; + + istart = 1; + jstart = 1; + kstart = 1; + + iend = cctk_lsh[0]-1; + jend = cctk_lsh[1]-1; + kend = cctk_lsh[2]-1; + + /* Do the evolution */ + + for (k=kstart; k<kend; k++) + { + for (j=jstart; j<jend; j++) + { + for (i=istart; i<iend; i++) + { + tmp[CCTK_GFINDEX3D(cctkGH,i,j,k)] = + 2*(1 - (dt2)*(1/dx2 + 1/dy2 + 1/dz2))* + phi[CCTK_GFINDEX3D(cctkGH,i,j,k)] - + phi_old[CCTK_GFINDEX3D(cctkGH,i,j,k)] + + (dt2) * + ( ( phi[CCTK_GFINDEX3D(cctkGH,i+1,j,k)] + +phi[CCTK_GFINDEX3D(cctkGH,i-1,j,k)] )/dx2 + +( phi[CCTK_GFINDEX3D(cctkGH,i,j+1,k)] + +phi[CCTK_GFINDEX3D(cctkGH,i,j-1,k)] )/dy2 + +( phi[CCTK_GFINDEX3D(cctkGH,i,j,k+1)] + +phi[CCTK_GFINDEX3D(cctkGH,i,j,k-1)] )/dz2); + } + } + } + + /* Update timeslices */ + for (k=kstart; k<kend; k++) + { + for (j=jstart; j<jend; j++) + { + for (i=istart; i<iend; i++) + { + phi_old[CCTK_GFINDEX3D(cctkGH,i,j,k)] + = phi[CCTK_GFINDEX3D(cctkGH,i,j,k)]; + phi[CCTK_GFINDEX3D(cctkGH,i,j,k)] + = tmp[CCTK_GFINDEX3D(cctkGH,i,j,k)]; + } + } + } + + /* Apply boundary conditions */ + WaveToyC_Boundaries(CCTK_PASS_CTOC); + + /* Synchronize */ + CCTK_SyncGroup(cctkGH,"wavetoy::scalarevolve"); + +} + + + + + + + + + |