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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
/*@@
@file SourceData.F77
@date
@author Gabrielle Allen
@desc
Elliptic initial data for wave equation
@enddesc
@@*/
#include "cctk.h"
#include "cctk_arguments.h"
#include "cctk_parameters.h"
subroutine UniformCharge(CCTK_FARGUMENTS)
c Find static field for a uniformly charge sphere
c
c That is, solve Nabla^2 phi = - 4 pi rho
c where rho = Q/(4/3 * Pi * R^3)
c where Q is the total charge and R is the sphere radius
implicit none
DECLARE_CCTK_FARGUMENTS
DECLARE_CCTK_PARAMETERS
CCTK_REAL pi
CCTK_REAL AbsTol(3), RelTol(3)
CCTK_REAL charge_factor
integer iphi,iMcoeff,iNcoeff
integer i,j,k,ierr
c Get variable indices for all grid functions
call CCTK_VarIndex (iMcoeff, "idscalarwaveelliptic::Mcoeff")
if (iMcoeff .lt. 0) then
call CCTK_WARN(0,"Grid variable index for Mcoeff not found")
end if
call CCTK_VarIndex (iNcoeff, "idscalarwaveelliptic::Ncoeff")
if (iNcoeff .lt. 0) then
call CCTK_WARN(0,"Grid variable index for Ncoeff not found")
end if
call CCTK_VarIndex (iphi,"wavetoy::phi")
if (iphi .lt. 0) then
call CCTK_WARN(0,"Grid variable index for iphi not found")
end if
c Set up all coefficients and initial guess for solution
pi = 4.0*atan(1.0)
charge_factor = 4.0d0*pi*charge*3.0d0/(4.0d0*pi*radius**3)
do k=1, cctk_lsh(3)
do j=1, cctk_lsh(2)
do i=1, cctk_lsh(1)
phi(i,j,k) = 0.0d0
Mcoeff(i,j,k) = 0.0d0
if (r(i,j,k) <= radius) then
Ncoeff(i,j,k) = charge_factor
else
Ncoeff(i,j,k) = 0d0
end if
end do
end do
end do
c Set tolerance for stopping elliptic solve
AbsTol(1)=1.0d-5
AbsTol(2)=1.0d-5
AbsTol(3)=1.0d-5
RelTol(1)=-1
RelTol(2)=-1
RelTol(3)=-1
c Set any options needed for the elliptic solve
c call Ell_SetStringKey(ierr, "yes", "EllLinFlat::Bnd::Robin")
c write (*,*) "SourceData: SetString ",ierr
call Ell_SetRealKey(ierr, 1.0d0, "EllLinFlat::Bnd::Robin::inf")
call Ell_SetIntKey (ierr, 1, "EllLinFlat::Bnd::Robin::falloff")
c Call elliptic solver to fill out phi
write (*,*) "charge: Going into elliptic solver"
call Ell_LinFlatSolver(
& ierr,
& cctkGH,
& iphi,
& iMcoeff, iNcoeff,
& AbsTol, RelTol,
& "sor")
write (*,*) "charge: Exit elliptic solver ierr =",ierr
c Set up last timestep ... assume time symmetry
do k=1, cctk_lsh(3)
do j=1, cctk_lsh(2)
do i=1, cctk_lsh(1)
phi_old(i,j,k) = phi(i,j,k)
end do
end do
end do
c Output exact solution if required
if (output_tmp==1) then
do k=1, cctk_lsh(3)
do j=1, cctk_lsh(2)
do i=1, cctk_lsh(1)
if (r(i,j,k) >= radius) then
temp(i,j,k) = charge/r(i,j,k)
else
temp(i,j,k) = charge/(2.0d0*radius**3)*
& (3.0d0*radius**2-r(i,j,k)**2)
end if
end do
end do
end do
call CCTK_OutputVarAsByMethod(ierr,cctkGH,
& "idscalarwaveelliptic::temp",
& "IOASCII_1D","phi_exact")
end if
return
end
|
