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
|
/*@@
@file WaveBinary.c
@date
@author Cactus Maintainers
@desc
Add a binary source term to the 3D scalar wave equation
@enddesc
@version $Header$
@@*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
static const char *rcsid = "$Header$";
CCTK_FILEVERSION(CactusWave_WaveBinarySource_WaveBinary_c)
/********************************************************************
********************* Local Data Types ***********************
********************************************************************/
/********************************************************************
********************* Local Routine Prototypes *********************
********************************************************************/
/********************************************************************
********************* Local Data *****************************
********************************************************************/
/********************************************************************
***************** External Routines Prototype ********************
********************************************************************/
void WaveBinaryC(CCTK_ARGUMENTS);
/********************************************************************
********************** External Routines ************************
********************************************************************/
void WaveBinaryC(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
static int firstcall=1;
const CCTK_REAL binary_pi=3.141592653589793;
CCTK_REAL charge_factor;
CCTK_REAL xsm,ysm,zsm,rad2m;
CCTK_REAL xsp,ysp,zsp,rad2p;
int i,j,k,vindex;
charge_factor = 3.0*binary_charge/(4.0*binary_pi*pow(binary_size,3));
/* Calculate the centers of the binary sources */
xsm = - binary_radius * cos(binary_omega*cctk_time);
ysm = - binary_radius * sin(binary_omega*cctk_time);
zsm = 0;
xsp = + binary_radius * cos(binary_omega*cctk_time);
ysp = + binary_radius * sin(binary_omega*cctk_time);
zsp = 0;
if ((CCTK_EQUALS(binary_verbose, "yes") && firstcall) ||
CCTK_EQUALS(binary_verbose, "debug"))
{
CCTK_VInfo(CCTK_THORNSTRING,
"Charge factor: %g\n",(double)charge_factor);
CCTK_VInfo(CCTK_THORNSTRING,
"Charge center (-): %g %g %g\n",
(double)xsm,(double)ysm,(double)zsm);
CCTK_VInfo(CCTK_THORNSTRING,
"Charge center (+): %g %g %g\n",
(double)xsp,(double)ysp,(double)zsp);
CCTK_VInfo(CCTK_THORNSTRING,
"Charge extent: %g\n",(double)binary_size);
}
firstcall=0;
for (k=0;k<cctk_lsh[2];k++)
{
for (j=0;j<=cctk_lsh[1];j++)
{
for (i=0;i<=cctk_lsh[0];i++)
{
vindex = CCTK_GFINDEX3D(cctkGH,i,j,k);
rad2m =
pow(x[vindex]-xsm,2) + pow(y[vindex]-ysm,2) + pow(z[vindex]-zsm,2);
rad2p =
pow(x[vindex]-xsp,2) + pow(y[vindex]-ysp,2) + pow(z[vindex]-zsp,2);
/* Note that both sources are positive, leading to a net
monopole moment */
if (rad2m<pow(binary_size,2))
{
phi[vindex] += 0.5*pow(CCTK_DELTA_TIME,2)*charge_factor;
}
if (rad2p<pow(binary_size,2))
{
phi[vindex] += 0.5*pow(CCTK_DELTA_TIME,2)*charge_factor;
}
}
}
}
/* Note that we do not need to sync anything, since each grid patch
has filled out its ghostzones */
}
|
