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
|
/*@@
@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"
#include "CactusElliptic/EllBase/src/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
CCTK_REAL pi
CCTK_REAL AbsTol(3), RelTol(3)
CCTK_REAL charge_factor
integer iphi,iMcoeff,iNcoeff
integer i,j,k,ierr
integer length
character*30 fsolver
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_SetStringKey(ierr, "yes", "EllLinFlat::Bnd::Robin")
call Ell_SetRealKey(ierr, 0d0, "EllLinFlat::Bnd::Robin::inf")
call Ell_SetIntKey (ierr, 1, "EllLinFlat::Bnd::Robin::falloff")
c Call elliptic solver to fill out phi
c Fill a fortran string with the name of the solver
call CCTK_FortranString(length,solver,fsolver)
call CCTK_INFO("Going into elliptic solver")
print *,"Solver:",fsolver(1:length),"."
call Ell_LinFlatSolver(
& ierr,
& cctkGH,
& iphi,
& iMcoeff, iNcoeff,
& AbsTol, RelTol,
& fsolver(1:length))
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
call CCTK_INFO("Leaving elliptic solver: solve failed")
end if
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_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
|