1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
/*@@
@file WaveToy.c
@date
@author Tom Goodale
@desc
Evolution routines for the wave equation solver
@enddesc
@version $Id$
@@*/
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
static const char *rcsid = "$Header$";
CCTK_FILEVERSION(CactusWave_WaveToyC_WaveToy_c);
void WaveToyC_Boundaries(CCTK_ARGUMENTS);
void WaveToyC_Evolution(CCTK_ARGUMENTS);
/*@@
@routine WaveToyC_Evolution
@date
@author Tom Goodale
@desc
Evolution for the wave equation
@enddesc
@calls CCTK_SyncGroup, WaveToyC_Boundaries
@@*/
void WaveToyC_Evolution(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
int i,j,k;
int vindex;
int istart, jstart, kstart, iend, jend, kend;
CCTK_REAL dx,dy,dz,dt,dx2,dy2,dz2,dt2;
CCTK_REAL dx2i,dy2i,dz2i;
CCTK_REAL factor;
/* 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;
dx2i = 1.0/dx2;
dy2i = 1.0/dy2;
dz2i = 1.0/dz2;
istart = 1;
jstart = 1;
kstart = 1;
iend = cctk_lsh[0]-1;
jend = cctk_lsh[1]-1;
kend = cctk_lsh[2]-1;
/* Do the evolution */
factor = 2*(1 - (dt2)*(dx2i + dy2i + dz2i));
for (k=kstart; k<kend; k++)
{
for (j=jstart; j<jend; j++)
{
for (i=istart; i<iend; i++)
{
vindex = CCTK_GFINDEX3D(cctkGH,i,j,k);
phi[vindex] = factor*
phi_p[vindex] - phi_p_p[vindex]
+ (dt2) *
( ( phi_p[CCTK_GFINDEX3D(cctkGH,i+1,j ,k )]
+phi_p[CCTK_GFINDEX3D(cctkGH,i-1,j ,k )] )*dx2i
+( phi_p[CCTK_GFINDEX3D(cctkGH,i ,j+1,k )]
+phi_p[CCTK_GFINDEX3D(cctkGH,i ,j-1,k )] )*dy2i
+( phi_p[CCTK_GFINDEX3D(cctkGH,i ,j ,k+1)]
+phi_p[CCTK_GFINDEX3D(cctkGH,i ,j, k-1)] )*dz2i);
}
}
}
}
/*@@
@routine WaveToyC_Boundaries
@date
@author Tom Goodale
@desc
Boundary conditions for the wave equation
@enddesc
@@*/
void WaveToyC_Boundaries(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
const char *bctype;
bctype = NULL;
if (CCTK_EQUALS(bound,"flat") || CCTK_EQUALS(bound,"static") ||
CCTK_EQUALS(bound,"radiation") || CCTK_EQUALS(bound,"robin") ||
CCTK_EQUALS(bound,"none"))
{
bctype = bound;
}
else if (CCTK_EQUALS(bound,"zero"))
{
bctype = "scalar";
}
/* Uses all default arguments, so invalid table handle -1 can be passed */
if (bctype && Boundary_SelectVarForBC (cctkGH, CCTK_ALL_FACES, 1, -1,
"wavetoy::phi", bctype) < 0)
{
CCTK_WARN (0, "WaveToyC_Boundaries: Error selecting boundary condition");
}
}
|
