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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
|
/*@@
@file SourceData.F77
@date
@author Gabrielle Allen
@desc
Elliptic initial data for wave equation
@enddesc
@version $Header$
@@*/
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"
#include "EllBase.h"
subroutine UniformCharge(CCTK_ARGUMENTS)
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_ARGUMENTS
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_FUNCTIONS
CCTK_REAL pi
CCTK_REAL AbsTol(3), RelTol(3)
CCTK_REAL charge_factor
integer iphi,iMcoeff,iNcoeff
integer i,j,k,ierr
CCTK_INT length
character*30 fsolver
character*200 infoline
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
call Ell_SetStrKey(ierr, "yes","EllLinFlat::Bnd::Robin")
call Ell_SetRealKey(ierr, 0d0, "EllLinFlat::Bnd::Robin::inf")
call Ell_SetIntKey (ierr, 1, "EllLinFlat::Bnd::Robin::falloff")
c Set parameters for specific solvers
if (CCTK_EQUALS(solver,"sor")) then
call Ell_SetIntKey(ierr, sor_maxit,"Ell::SORmaxit")
end if
call CCTK_FortranString(length,solver,fsolver)
write(infoline,'("Going into elliptic solver ",A)') fsolver
call CCTK_INFO(infoline)
c Call elliptic solver to fill out phi
call Ell_LinFlatSolver(
& ierr,
& cctkGH,
& iphi,
& iMcoeff, iNcoeff,
& AbsTol, RelTol,
& fsolver)
if (ierr .eq. ELL_SUCCESS) then
call CCTK_INFO("Leaving elliptic solver: solve successful")
else if (ierr .eq. ELL_NOCONVERGENCE) then
call CCTK_INFO("Leaving elliptic solver: solver failed to converge")
else if (ierr .eq. ELL_NOSOLVER) then
call CCTK_INFO("Elliptic solver not found")
else
write(infoline,'("Leaving elliptic solver: solve failed (Error ",I3,")")') ierr
call CCTK_INFO(infoline)
end if
c Set up last timestep ... assume (first order) time symmetry
do k=1, cctk_lsh(3)
do j=1, cctk_lsh(2)
do i=1, cctk_lsh(1)
phi_p(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
|
