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
|
/*@@
@file SourceData.F90
@date
@author Gabrielle Allen
@desc
Elliptic initial data for wave equation
Originally written in F77 fixed format. Converted to F90
free format by Peter Diener.
@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)
! Find static field for a uniformly charge sphere
!
! That is, solve Nabla^2 phi = - 4 pi rho
! where rho = Q/(4/3 * Pi * R^3)
! 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
! 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
! Set up all coefficients and initial guess for solution
pi = 4.0d0*atan(1.0d0)
charge_factor = 4.0d0*pi*charge*3.0d0/(4.0d0*pi*radius**3)
phi = 0.0d0
Mcoeff = 0.0d0
where ( r <= radius )
Ncoeff = charge_factor
elsewhere
Ncoeff = 0.0d0
end where
! 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
! Set any options needed for the elliptic solve
call Ell_SetStrKey (ierr, "yes", "EllLinFlat::Bnd::Robin")
call Ell_SetRealKey(ierr, 0.0d0, "EllLinFlat::Bnd::Robin::inf")
call Ell_SetIntKey (ierr, 1, "EllLinFlat::Bnd::Robin::falloff")
! 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)
! 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
! Set up last timestep ... assume (first order) time symmetry
phi_p = phi
! Output exact solution if required
if (output_tmp .eq. 1) then
where ( r .ge. radius )
temp = charge / r
elsewhere
temp = charge/(2.0d0*radius**3)*(3.0d0*radius**2-r**2)
end where
call CCTK_OutputVarAsByMethod(ierr,cctkGH, &
"idscalarwaveelliptic::temp", &
"IOASCII_1D","phi_exact")
end if
return
end
|
