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
|
/*@@
@file WaveToy.cc
@date
@author Tom Goodale
@desc
Evolution routines for the wave equation solver
@enddesc
@version $Id$
@@*/
#include <stddef.h>
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#define val(gridfunc,i,j,k) gridfunc[CCTK_GFINDEX3D(cctkGH,i,j,k)]
/*@@
@routine WaveToyC_Evolution
@date
@author Tom Goodale
@desc
Evolution for the wave equation
@enddesc
@@*/
extern "C" void WaveToyCXX_Evolution(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
// Set up shorthands
CCTK_REAL dx = CCTK_DELTA_SPACE(0);
CCTK_REAL dy = CCTK_DELTA_SPACE(1);
CCTK_REAL dz = CCTK_DELTA_SPACE(2);
CCTK_REAL dt = CCTK_DELTA_TIME;
CCTK_REAL dx2 = dx*dx;
CCTK_REAL dy2 = dy*dy;
CCTK_REAL dz2 = dz*dz;
CCTK_REAL dt2 = dt*dt;
CCTK_REAL dx2i = 1.0/dx2;
CCTK_REAL dy2i = 1.0/dy2;
CCTK_REAL dz2i = 1.0/dz2;
CCTK_REAL factor = 2*(1 - (dt2)*(dx2i + dy2i + dz2i));
int istart = 1;
int jstart = 1;
int kstart = 1;
int iend = cctk_lsh[0]-1;
int jend = cctk_lsh[1]-1;
int kend = cctk_lsh[2]-1;
//
// Do the evolution
//
for (int k=kstart; k<kend; k++)
{
for (int j=jstart; j<jend; j++)
{
for (int i=istart; i<iend; i++)
{
int vindex = CCTK_GFINDEX3D(cctkGH,i,j,k);
phi[vindex] =
factor*phi_p[vindex] - phi_p_p[vindex]
+ dt2 *
( ( val( phi_p, i+1,j ,k ) + val( phi_p, i-1,j ,k) )*dx2i
+( val( phi_p, i ,j+1,k ) + val( phi_p, i ,j-1,k) )*dy2i
+( val( phi_p, i ,j ,k+1) + val( phi_p, i ,j, k-1) )*dz2i
);
}
}
}
}
/*@@
@routine WaveToyC_Boundaries
@date
@author Tom Goodale
@desc
Boundary conditions for the wave equation
@enddesc
@@*/
extern "C" void WaveToyCXX_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,"WaveToyCXX_Boundaries: Error selecting boundary condition");
}
}
|
