-
Notifications
You must be signed in to change notification settings - Fork 0
/
potsetup.F90
267 lines (228 loc) · 8.04 KB
/
potsetup.F90
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
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
subroutine potsetup
implicit none
include 'runhydro.h'
include 'pot.h'
!********************************************************************
!*
! potsetup does intitialization for the poisson solver
! package. In order of occurance this is what potsetup
! does:
!
! -> initialize trigonometric functions used by bessel
!
! -> call tm and sm to fill in the Green function tmr and smz
! respectively
!
! -> setup the geometric arrays for the ADI solver
!
!*
!********************************************************************
!*
!* Global Variables
real, dimension(numr,numz,numr,mmax) :: tmr
real, dimension(numr,numz,numz,mmax) :: smz
common /green_functions/ tmr, smz
real, dimension(numphi,mmax) :: bes_cos, bes_sin
common /bessel_trig/ bes_cos, bes_sin
real, dimension(numr) :: ar, cr, alphar
real, dimension(numr,numphi) :: brb
common /ADI_R_sweep/ ar, cr, alphar, brb
real :: az, cz
real, dimension(numr) :: alphaz, betaz
real, dimension(numz) :: bzb
real, dimension(numr,numphi) :: elambdazb
common /ADI_Z_sweep/ az, cz, alphaz, betaz, bzb, elambdazb
real :: gamma, piinv, four_pi
common /pot_constants/ gamma, piinv, four_pi
real :: grav
common /constants/ grav
real :: dr, dz, dphi, drinv, dzinv, dphiinv
common /coord_differentials/ dr, dz, dphi, drinv, dzinv, dphiinv
real, dimension(numr) :: rhf, r, rhfinv, rinv
real, dimension(numz) :: zhf
real, dimension(numphi) :: phi
common /grid/ rhf, r, rhfinv, rinv, zhf, phi
real, dimension(numr) :: rhf_g, r_g, rhfinv_g, rinv_g
real, dimension(numz) :: zhf_g
common /global_grid/ rhf_g, r_g, rhfinv_g, rinv_g, zhf_g
real, dimension(numr,numz,numphi) :: pot, rho
common /poisson/ pot, rho
! integer, dimension(3) :: boundary_condition
! common /boundary_conditions/ boundary_condition
!*
!********************************************************************
!*
!* Local Variables
real :: mindex, lindex, drinv2, dphiinv2
real, dimension(numphi) :: elm, m1mode
real, dimension(numr) :: betar
integer :: j, k, l, m, lstop, index, mode
!*
!********************************************************************
! initialize the local variables
mindex = 0.0
lindex = 0.0
drinv2 = 0.0
dphiinv2 = 0.0
elm = 0.0
m1mode = 0.0
betar = 0.0
j = 0
k = 0
l = 0
m = 0
lstop = 0
index = 0
mode = 0
piinv = 1.0 / pi
four_pi = 4.0 * pi
! write(*,*) 'POTSETUP: pi ',pi
! write(*,*) 'POTSETUP: piinv ',piinv
! write(*,*) 'POTSETUP: four_pi ', four_pi
! setup the trig arrays used by bessel, use lindex and mindex to
! avoid casts which are relatively expensive
mindex = 1.0
lindex = 1.0
if( isym == 3 ) then
do m = 1, mmax
do l = philwb, phiupb
bes_cos(l,m) = cos(dphi*(mindex-1.0)*(2.0*lindex-1.0))
bes_sin(l,m) = sin(dphi*(mindex-1.0)*(2.0*lindex-1.0))
lindex = lindex + 1.0
enddo
lindex = 1.0
mindex = mindex + 1.0
enddo
else
do m = 1, mmax
do l = philwb, phiupb
bes_cos(l,m) = cos(0.5*dphi*(mindex-1.0)*(2.0*lindex-1.0))
bes_sin(l,m) = sin(0.5*dphi*(mindex-1.0)*(2.0*lindex-1.0))
lindex = lindex + 1.0
enddo
lindex = 1.0
mindex = mindex + 1.0
enddo
endif
! open(unit=10,file='bes_sin',form='unformatted',status='unknown')
! write(10) bes_sin
! close(10)
! open(unit=11,file='bes_cos',form='unformatted',status='unknown')
! write(11) bes_cos
! close(11)
! call tm and sm to fill in the Green functions for the top (and bottom)
! potential - tmr and the Green function for the side potential - smz
call tm(tmr)
call sm(smz)
! open(unit=10,file='/scratch1/phmotl/tmr',form='unformatted',status='unknown')
! read(10) tmr
! close(10)
! open(unit=11,file='/scratch1/phmotl/smz',form='unformatted',status='unknown')
! read(11) smz
! close(11)
! the rest of the code in potsetup sets up arrays that hold geometric
! information for the discretization of Poisson's equation. The
! arrays are all used by the ADI solver
gamma = dzinv * dzinv
drinv2 = drinv * drinv
dphiinv2 = dphiinv * dphiinv
! write(*,*) 'POTSETUP: gamma ',gamma
! write(*,*) 'POSTETUP: drinv2 ',drinv2
! write(*,*) 'POTSETUP: dphiinv2 ',dphiinv2
lstop = numphi_by_two + 1
do l = philwb, phiupb
if( isym == 3 ) then
if( l <= lstop ) then
mode = (l-1)*(l-1)
else
mode = (l-lstop)*(l-lstop)
endif
else
if( l <= lstop ) then
mode = l - 1
else
mode = l - lstop
endif
endif
m1mode(l) = (-1.0)**mode
elm(l) = cos(mode*dphi)
enddo
! open(unit=10,file='m1mode',form='unformatted',status='unknown')
! write(10) m1mode
! close(10)
! open(unit=11,file='elm',form='unformatted',status='unknown')
! write(11) elm
! close(11)
do j = rlwb, rupb
alphar(j) = -r_g(j+1)*rhfinv_g(j)*drinv2
betar(j) = -r_g(j)*rhfinv_g(j)*drinv2
enddo
! open(unit=12,file='alphar',form='unformatted',status='unknown')
! write(12) alphar
! close(12)
! open(unit=13,file='betar',form='unformatted',status='unknown')
! write(13) betar
! close(13)
! have to used indexed global radius array to initialize
! alphaz and betaz because they will be used when the radial
! data is block distributed across numz_procs and numr_dd
! is not in general equal to numr_dd_pad
do j = rlwb, rupb
alphaz(j) = -r_g(j+1)*rhfinv_g(j)*drinv2
betaz(j) = -r_g(j)*rhfinv_g(j)*drinv2
enddo
! open(unit=14,file='alphaz',form='unformatted',status='unknown')
! write(14) alphaz
! close(14)
! open(unit=15,file='betaz',form='unformatted',status='unknown')
! write(15) betaz
! close(15)
ar = betar
cr = alphar
az = - gamma
cz = - gamma
do l = philwb, phiupb
do j = rlwb + 1, rupb
brb(j,l) = 2.0*drinv2 - 2.0*(elm(l)-1.0)*dphiinv2*rhfinv_g(j)*rhfinv_g(j)
enddo
enddo
if( isym == 3 ) then
do l = philwb, phiupb
brb(2,l) = -alphar(2) - 2.0*betar(2) - 2.0*(elm(l)-1.0)*dphiinv2*rhfinv_g(2)*rhfinv_g(2)
enddo
else
do l = philwb, phiupb
brb(2,l) = -alphar(2) + (m1mode(l)-1.0)*betar(2) -2.0*(elm(l)-1.0)*dphiinv2*rhfinv_g(2)*rhfinv_g(2)
enddo
endif
! open(unit=17,file='brb',status='unknown',form='unformatted')
! write(17) brb
! close(17)
do k = zlwb, zupb
bzb(k) = 2.0*gamma
enddo
if( isym /= 1 ) bzb(2) = gamma
! open(unit=18,file='bzb',form='unformatted',status='unknown')
! write(18) bzb
! close(18)
do l = philwb, phiupb
do j = rlwb, rupb
elambdazb(j,l) = -2.0*drinv2 + 2.0*(elm(l)-1.0)*dphiinv2*rhfinv_g(j)*rhfinv_g(j)
enddo
enddo
if( isym == 3 ) then
do l = philwb, phiupb
elambdazb(2,l) = alphaz(2) + 2.0*betaz(2) + 2.0*(elm(l)-1.0)*dphiinv2*rhfinv(2)*rhfinv(2)
enddo
endif
if( isym /= 3 ) then
do l = philwb, phiupb
elambdazb(2,l) = alphaz(2) - (m1mode(l)-1.0)*betaz(2) +2.0*(elm(l)-1.0)*dphiinv2*rhfinv(2)*rhfinv(2)
enddo
endif
! open(unit=19,file='elambdazb',form='unformatted',status='unknown')
! write(19) elambdazb
! close(19)
! print*,"mmax potsetup",mmax
return
end subroutine potsetup