-
Notifications
You must be signed in to change notification settings - Fork 3
/
piccom.f
234 lines (221 loc) · 9.24 KB
/
piccom.f
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
integer npartmax,npart,nr,nth,npsi,ndim,np
c Number of particles: npartmax, radial and theta mesh size: nr, nth.
c Don't change anything else.
parameter (npartmax=400000,np=1,ndim=6)
c Use of particle advance subcycling in inner regions for accuracy.
logical lsubcycle
c Integrator type. True=old, False=new symplectic schemes
logical verlet
c CIC definitions
logical LCIC,collcic
integer NRUSED,NTHUSED,NPSIUSED,NRFULL,NTHFULL,NPSIFULL
parameter (LCIC=.true.)
integer nrsize,nthsize,npsisize
c These correspond to nrfull, nthfull.and npsifull
parameter (nrsize=306,nthsize=31,npsisize=31)
c Positions and velocities of particles (6-d phase-space).
real xp(ndim,npartmax)
real vzinit(npartmax)
real dtprec(npartmax)
c Flag of particle slot status (e.g. in use or not)
integer ipf(npartmax)
c The potential normalized to Te/e
real phi(0:nrsize,0:nthsize,0:npsisize)
c The potential on axis (cos(theta)=+-1) before averaging
real phiaxis(0:nrsize,2,0:npsisize)
c Charge density
real rho(0:nrsize,0:nthsize,0:npsisize)
real rhoDiag(0:nrsize,0:nthsize,0:npsisize)
c Injection complement. How many particles to reinject each step
integer ninjcomp
c Highest occupied particle slot.
integer iocprev
real pi
parameter (pi=3.1415927)
real cerr,bdyfc,Ti,vd,cd,cB,Bz
logical diags,lplot,ldist,linsulate,lfloat,lat0,lap0,localinj
logical lfixedn
integer myid,numprocs
real rmtoz
common /piccom/xp,npart,vzinit,dtprec,phi,rho,rhoDiag,cerr,bdyfc
$ ,Ti,vd,cd,cB,diags,ninjcomp,lplot,ldist,linsulate,lfloat
$ ,lat0,lap0 ,localinj,lfixedn,myid,numprocs,rmtoz,ipf,iocprev
$ ,Bz,lsubcycle,verlet,collcic,phiaxis
c *******************************************************************
c Momenta of the distribution function
c The particle number
real psum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
c The sum of particle radial velocities
real vrsum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
c Sum of theta velocities
real vtsum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
c Sum of phi velocities
real vpsum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
c The sum of particle velocities squared
real vr2sum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vt2sum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vp2sum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
c Off diagonal terms
real vrtsum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vrpsum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vtpsum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
c The sum of particle xyz-velocities
real vxsum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vysum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vzsum(1:nrsize-1,1:nthsize-1,1:npsisize-1)
c Diagnostic sums
real pDiag(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vrDiag(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vtDiag(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vpDiag(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vr2Diag(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vt2Diag(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vp2Diag(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vrtDiag(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vrpDiag(1:nrsize-1,1:nthsize-1,1:npsisize-1)
real vtpDiag(1:nrsize-1,1:nthsize-1,1:npsisize-1)
c Total sum of particle xyz-velocities, i.e. total current curr(4)
c is the particle sum
real curr(4)
common /momcom/psum,vrsum,vtsum,vpsum,vr2sum,vt2sum,vp2sum ,vrtsum
$ ,vrpsum,vtpsum,vzsum,vxsum,vysum,curr,pDiag,vrDiag,vtDiag
$ ,vpDiag,vr2Diag,vt2Diag,vp2Diag ,vrtDiag,vrpDiag,vtpDiag
c*********************************************************************
c Radius mesh
real r(0:nrsize),rcc(0:nrsize)
c Theta angle mesh
real th(0:nthsize),tcc(0:nthsize)
c Theta mesh radians
real thang(0:nthsize)
c Poloidal (Psi) Mesh. Only cell center are useful
real pcc(0:npsisize)
c Inverse of volume of radial shells.
real volinv(0:nrsize)
c Precalculation functions
integer nrpre,ntpre,nppre
parameter (nrpre=4*nrsize,ntpre=4*nthsize,nppre=4*npsisize)
integer irpre(nrpre),itpre(ntpre),ippre(nppre)
real rfac,tfac,pfac
c Non-uniform handling quantities.
real hr(0:nrsize+1),zeta(0:nrsize+1),zetahalf(0:nrsize+1),
$ cminus(nrsize),cmid(nrsize),cplus(nrsize)
c Lower limit of averaging range. 0.6 by default
real avelim
c Parallel or serial solving
logical cgparallel
c Parallel bloc solver arguments
integer idim1,idim2,idim3
common /meshcom/r,rcc,th,tcc,thang,volinv,irpre,itpre,rfac,tfac,
$ pcc,ippre,pfac, hr,zeta,zetahalf,cminus,cmid,cplus ,avelim
$ ,nr,NRFULL,NRUSED,NPSIFULL,NPSIUSED,nth,npsi,NTHFULL,NTHUSED
$ ,cgparallel,idim1,idim2,idim3
c********************************************************************
c Random interpolate data.
integer nvel,nQth
parameter (nvel=50)
parameter (nQth=200)
real Qcom(nQth)
real Gcom(nvel,nQth)
real Vcom(nvel)
real pu1(nvel),pu2(nvel)
real Pc(nQth,nvel)
c New BC
integer bcphi,bcr
logical infdbl
c Reinjection flux as a function of cos(theta) (line) and chi (column,
c from 0 to 9)
common /rancom/Gcom,Vcom,Qcom,pu1,pu2,Pc,infdbl,bcphi,bcr
c********************************************************************
c diagnostic data
integer nvmax,nrein,nreintry,ninner,nstepmax
parameter (nvmax=60,nstepmax=10001)
real nvdiag(nvmax),nvdiagave(nvmax),vdiag(nvmax)
real vrdiagin(nvmax),vtdiagin(nvmax)
real vrange
real diagrho(nrsize),diagphi(nrsize)
real diagchi(0:nthsize)
real phiout
integer partz,fieldz,epressz,enccharge,lorentz
parameter(enccharge=1,fieldz=2,epressz=3,partz=4,lorentz=5)
c Total particle flux to probe
real fluxprobe(nstepmax)
c Total momentum flux to probe
real zmomprobe,xmomprobe,ymomprobe
c Total energy collected
real enerprobe
c Z-momentum flux across outer boundary.
real zmout,xmout,ymout
c Combined zmom data: fields, electron pressure, ion momentum.
c For inner 1, outer 2. zmom also carries the probe charge
real zmom(nstepmax,5,2),xmom(nstepmax,2:5,2),ymom(nstepmax,2:5,2)
c enertot is the reduced enerprobe for each time-step
real enertot(nstepmax)
c Number of particles striking probe in theta/psi cell
real nincellstep(nthsize,npsisize,0:nstepmax)
real vrincellstep(nthsize,npsisize,0:nstepmax)
real vr2incellstep(nthsize,npsisize,0:nstepmax)
real nincell(nthsize,npsisize)
real vrincell(nthsize,npsisize)
real vr2incell(nthsize,npsisize)
c Impose the bohm condition or not
logical bohm
c Ave flux
real fincellave(nthsize,npsisize)
c Ave radial mom flux
real vrincellave(nthsize,npsisize)
c Ave radial vr2 flux
real vr2incellave(nthsize,npsisize)
c Number of particles reinjected per theta cell.
c integer noutrein(nth),ivoutrein(nth)
c Sum and average of potentials at which particles were reinjected.
real spotrein,averein
c Flux of reinjection.
real fluxrein
c Density at infinity.
real rhoinf
c Number of trapped reinjections
integer ntrapre
c Coefficient of density deficit, for external solution
real adeficit
c Cell in which to accumulate distribution functions
integer ircell,itcell
common /diagcom/nvdiag,nvdiagave,vdiag,vrange,diagrho,diagphi,
$ diagchi,phiout,nrein,nreintry,ninner,fluxprobe,nincellstep
$ ,vrincellstep,vr2incellstep,nincell,vrincell,vr2incell,rhoinf
$ ,vrdiagin,vtdiagin, spotrein,averein ,fluxrein ,ntrapre
$ ,adeficit, ircell,itcell ,zmout,xmout,ymout ,zmomprobe
$ ,ymomprobe,xmomprobe,fincellave ,vrincellave,vr2incellave
$ ,zmom,xmom,ymom ,enerprobe ,enertot, bohm
c*********************************************************************
c Poisson coefficients for iterative solution, etc.
real debyelen,vprobe,Ezext
real apc(0:nrsize),bpc(0:nrsize)
real cpc(0:nrsize,0:nthsize),dpc(0:nrsize,0:nthsize)
real epc(0:nrsize,0:nthsize)
real fpc(0:nrsize,0:nthsize)
real gpc(0:nthsize,0:npsisize,1:5)
c Flag indicating to use biconjugate gradient method (not min. res.)
logical lbcg
common /poisson/debyelen,vprobe,Ezext,apc,bpc,cpc,dpc,fpc,epc,gpc
$ ,lbcg
c*********************************************************************
c Smoothing steps
integer nstepsave,nsamax,diagsamp
logical samp
common /stepave/nstepsave,nsamax,diagsamp,samp
c*********************************************************************
c Orbit plotting storage for tracking the first norbits orbits.
integer nobsmax,norbits
parameter (nobsmax=100)
real xorbit(nstepmax,nobsmax),yorbit(nstepmax,nobsmax),
$ zorbit(nstepmax,nobsmax)
real vxorbit(nstepmax,nobsmax),vyorbit(nstepmax,nobsmax),
$ vzorbit(nstepmax,nobsmax),rorbit(nstepmax,nobsmax)
integer iorbitlen(nobsmax)
common /orbits/norbits,iorbitlen,xorbit,yorbit,zorbit,rorbit,
$ vxorbit,vyorbit,vzorbit
c*********************************************************************
c Data necessary for the orbit tracking
logical orbinit
integer maxsteps,trackinit
common /orbtrack/orbinit,maxsteps,trackinit