PPT - University of California, Irvine
Download
Report
Transcript PPT - University of California, Irvine
GTC – Gyrokinetic Toroidal Code
Training course: Part I
Ihor Holod
University of California, Irvine
Outline
• Introduction
• Gyrokinetic approach
• GK equations of motion. GK Poisson’s equation. Delta-f method
• Fluid-kinetic hybrid electron model: continuity equation,
Ampere’s law, drift-kinetic equation
• Tokamak geometry. Magnetic coordinates
• Particle-in-cell approach. General structure of GTC. Grid.
Parallel computing.
• Setup and loading
March 5-7, 2011
GTC Training Course - PKU
2
Introduction
• Gryrokinetic simulations – important tool to study turbulence
in magnetically confined plasma
• Main target – low frequency microinstabilities
• Source of free energy is needed to drive instabilities in plasma
• Pressure gradient of ions, electrons or fast (energetic) particles
• Ion temperature gradient mode (ITG)
• Trapped electron mode (TEM)
• Electron-temperature gradient mode (ETG)
• Energetic particle mode (EPM)
March 5-7, 2011
GTC Training Course - PKU
3
Gyrokinetic approach
• Main advantage: gyroaveraging reduces dimensionality of
particle dynamics from 6 to 5 [Hahm et al. PoF’88]
• Gyrokinetic ordering:
• Frequency smaller than gyrofrequency
• Small fluctuation amplitude
e B f
~
~
~
~ 1
T
B
F
• Additionally
k s ~ 1
k||
k
~ 1
• Magnetic field inhomogeneity is small:
March 5-7, 2011
GTC Training Course - PKU
ln B ~ B 1
4
Charged particle motion in constant uniform
magnetic filed
March 5-7, 2011
GTC Training Course - PKU
5
Guiding center motion in the inhomogeneous
magnetic filed
dv||
dX
ˆ
ˆ
v||b 0 v d
b 0 B0
dt
dt
m
• Magnetic drift velocity
v d vc v g
d
0
dt
• Magnetic curvature drift
vc
v||2
bˆ 0
• Magnetic gradient drift
vg
ˆ
b 0 B0
m
ing
x,v XGC , , v|| , Gyroaverag
XGY , , v||
March 5-7, 2011
GTC Training Course - PKU
6
Gyrocenter equations of motion
• Gyrocenter equations of motion [Brizard and Hahm, RMP‘07]:
dX
v||bˆ 0 v d v E
dt
1 B*
Z A||
B0 Z
dt
m B0
mc t
dv||
• Modified magnetic field:
B B B B 0
*
*
0
B0v||
bˆ 0 B
• Magnetic field perturbation:
B B A||bˆ 0
• ExB drift velocity:
cbˆ 0
vE
B0
March 5-7, 2011
GTC Training Course - PKU
7
Gyrokinetic Vlasov-Maxwell’s system of equations
• Gyrokinetic Vlasov equation
f
f
X f v||
0
t
v||
• Gyrokinetic Poisson’s equation
4Zi2 ni
~
4 Zi ni ene
Ti
• Gyrokinetic Ampere’s law
2 A||
4
eneu||e Zi niu||i
c
• Charge density and current
Zn Z
B0
m
Znu|| Z
March 5-7, 2011
dv d f
||
B0
m
dv d v f
||
||
GTC Training Course - PKU
8
Particle-in-cell method
• Phase space is randomly sampled by Lagrangian “markers”
• Number of markers is much smaller than # of physical particles
• Dynamics of markers is governed by the same equations of motion as for
physical particles
• Fluid moments, like density and current, are calculated on Eulerian grid
Monte-Carlo sampling of phase
space volume
March 5-7, 2011
Fluid moments and potentials
are calculated on stationary grid
GTC Training Course - PKU
9
Delta-f method
• Discrete particle noise – statistical error from MC sampling
• Noise can be reduced by increasing # of markers, hence increasing CPU time
• Delta-f method to reduce noise:
Perturbed part of the distribution function is considered
as additional dynamical quantity associated with particle,
and is governed by gyrokinetic equation
Fluid moments are calculated integrating contributions from
each individual particle
dX
v||bˆ 0 v d v E
dt
1 B*
Z A||
B0 Z
dt
m B0
mc t
dv||
B
B
A|| 1 ln f 0
B*
dw
1
1 w v||
v E ln f 0
B0 Z
m v
dt
B
B
B
c
t
0
0
||
0
March 5-7, 2011
GTC Training Course - PKU
10
Fluid-kinetic hybrid electron model - I
• Simulations of various Alfven modes
• Electron continuity equation
n u
ne
n
B0bˆ 0 0 ||e B0 v E 0 n0 v* v E ln B0 0
t
B0
B0
• Perturbed diamagnetic drift velocity
v*
1
bˆ 0 P|| P
n0 me e
• Perturbed pressure
B0
dv|| d B0f
m
B
P|| 0 dv|| d m v||2f
m
P
March 5-7, 2011
GTC Training Course - PKU
11
Fluid-kinetic hybrid electron model - II
• Lowest order solution of electron DKE based on /k||v|| expansion
f e(0)
f0
( 0)
eeff
Te
ln f 0
ln n0
• Electron adiabatic response
( 0)
eeff
Te
ne
n0
• Inverse gyrokinetic Ampere’s law
4
4
en eu||e 2 A||
Z i ni u||i
c
c
• Faraday’s law (linear tearing mode is removed)
1 A||
E|| bˆ 0 bˆ 0 eff
c t
March 5-7, 2011
GTC Training Course - PKU
12
Fluid-kinetic hybrid electron model - III
• High-order electron kinetic response
dwe f e( 0 )
f e( 0 )
1
we v E ln f 0 e
dt
f0
t
f
0
f e( 0 ) e
f e( 0) ev|| A||
v d
v E
Te
f0
cTe t
f0
(1)
eeff
Te
ne(1)
n0
• Zonal flow – nonlinearly self-generated large scale ExB flow with
electrostatic potential having k||=0
March 5-7, 2011
GTC Training Course - PKU
13
Tokamak geometry
• Equilibrium magnetic field (covariant representation) in flux coordinates
B0 g I
• Helicity of magnetic field is determined by safety factor q (number of
toroidal turns of magnetic field line per one poloidal turn)
d
1
d q
Magnetic flux surface
March 5-7, 2011
GTC Training Course - PKU
14
Magnetic coordinates - definition
• Covariant representation of the magnetic field
B0 g I
• Contravariant representation of the magnetic field
B0 q
• Jacobian
J
1
B02 ( , )
g ( )q( ) I ( )
B0 1 cos O 2
g 1 O
I 0O 2
2
0 O
0 O
0 O
March 5-7, 2011
GTC Training Course - PKU
15
Magnetic coordinates - example
• Parallel derivative
B 0 q
q 1
q
J q
J 1
March 5-7, 2011
GTC Training Course - PKU
16
GTC grids
• Two spatial grids
• Coarse grid for equilibrium and profiles
• Field-aligned simulation mesh for fields
• Spline interpolation is used to calculate equilibrium
quantities on simulation grid and particle positions
March 5-7, 2011
GTC Training Course - PKU
17
Simulation grid
• Radial grid in (i=0:mpsi)
• Poloidal angle () grid with continuous indexing
(ij=1:mgrid) over the whole poloidal plane
• Number of poloidal grid points mtheta(i)differs
from one flux surface to another
ij
igrid(0)=1
do i=1,mpsi
igrid(i)=igrid(i-1) + mtheta(i-1) + 1
enddo
March 5-7, 2011
2
mtheta(i)
j integer
GTC Training Course - PKU
ij igrid(i) j
18
Parallel computing
• Distribution of independent tasks between different processors
• Domain decomposition: each processor simulates it’s own part of
simulation domain
• Particle decomposition: each processor is responsible for it’s own
set of particles
Supercomputer structure
March 5-7, 2011
Core 0 Core 1
Core 0 Core 1
Core 2 Core 3
Core 2 Core 3
Node 0
Node 1
GTC Training Course - PKU
19
MPI-Parallelization
• Parallelization between nodes
• Each node has independent memory
• Inter-node communication through special commands
MPI_BCAST(send_buff,buff_size,data_type,send_pe,comm,ierror)
MPI_REDUCE(send_buff,recv_buff,buff_size,data_type,operation,
recv_pe,comm,ierror)
MPI_ALLREDUCE(send_buff,recv_buff,buff_size,data_type,operati
on,comm,ierror)
March 5-7, 2011
GTC Training Course - PKU
20
OpenMP-Parallelization
• Parallelization between cores of one node
• Each core has access to common memory
• Used for parallelization of cycles
• Useful with modern multi-core processors (GPU)
!$omp parallel do private(m)
do m=1,mi
…
enddo
March 5-7, 2011
GTC Training Course - PKU
21
GTC structure
Dynamics
ne
ge1&fi
A||
Fields
ue
A|| ZF es
ind
Sources
A|| ui
March 5-7, 2011
ne1 ni ue1 ne
GTC Training Course - PKU
22
main.F90
CALL INITIAL(cdate,timecpu)
CALL SETUP
CALL LOAD
CALL LOCATEI
if(fload>0)CALL LOCATEF
! main time loop
do istep=1,mstep
do irk=1,2
! idiag=0: do time history diagnosis at irk=1 & istep=ndiag*N
idiag=mod(irk+1,2)+mod(istep,ndiag)
if(idiag==0)then
if(ndata3d==1)CALL DATAOUT3D !write 3D fluid data
if(track_particles==1)CALL PARTOUT !write particle data
endif
CALL FIELD_GRADIENT
if(magnetic==1)CALL PUSHFIELD
CALL PUSHI
CALL SHIFTI(mimax,mi)
CALL LOCATEI
CALL CHARGEI
if(fload>0)then
CALL PUSHF
CALL SHIFTF(mfmax,mf)
CALL LOCATEF
CALL CHARGEF
endif
if(magnetic==0)then
CALL POISSON_SOLVER("adiabatic-electron")
else
CALL POISSON_SOLVER("electromagnetic")
CALL FIELD_SOLVER
endif
March 5-7, 2011
GTC Training Course - PKU
23
main.F90 (cont)
! push electron, sub-cycling
do i=1,ncycle*irk
! 1st RK step
CALL PUSHE(i,1)
CALL SHIFTE(memax,me)
! 2nd RK step
CALL PUSHE(i,2)
CALL SHIFTE(memax,me)
! electron-ion and electron-electron collision
if(ecoll>0 .and. mod(i,ecoll)==0)then
call lorentz
call collisionee
endif
enddo
CALL CHARGEE
CALL POISSON_SOLVER("kinetic-electron")
enddo
if(idiag==0)CALL DIAGNOSIS
enddo
! i-i collisions
if(icoll>0 .and. mod(istep,icoll)==0)call collisionii
if(mod(istep,mstep/msnap)==0 .or. istep*(1-irun)==ndiag)CALL SNAPSHOT
if(mod(istep,mstep/msnap)==0)CALL RESTART_IO("write")
enddo
CALL FINAL(cdate,timecpu)
March 5-7, 2011
GTC Training Course - PKU
24
analytic.F90
• Specifies profiles for different analytic equilibriums
if(numereq==1)then ! Cyclone case
psiw=0.0375 !poloidal flux at wall
! q and zeff profile is parabolic: q=q1+q2*psi/psiw+q3*(psi/psiw)^2
q=(/0.82,1.1,1.0/) !cyclone case
er=(/0.0,0.0,0.0/)
itemp0=1.0
! on-axis thermal ion temperature, unit=T_e0
ftemp0=2.0
! on-axis fast ion temperature, unit=T_e0
fden0=1.0e-5
! on-axis fast ion density, unit=n_e0
! density and temperature profiles are hyperbolic:
! ne=1.0+ne1*(tanh(ne2-(psi/psiw)/ne3)-1.0)
ne=(/0.205,0.30,0.4/) !Cyclone base case
te=(/0.415,0.18,0.4/)
ti=(/0.415,0.18,0.4/) ! ITG
tf=ti ! fast ion temperature profile
nf=ne ! fast ion density profile
March 5-7, 2011
GTC Training Course - PKU
25
eqdata.F90
• Equilibrium and profiles are represented using 2D-spline on (psi,
theta) or 1D-spline on (psi)
! allocate memory for equilibrium spline array
allocate (stpp(lsp),mipp(lsp),mapp(lsp),nepp(3,lsp),tepp(3,lsp),&
nipp(3,lsp),tipp(3,lsp),zepp(3,lsp),ropp(3,lsp),erpp(3,lsp),&
qpsi(3,lsp),gpsi(3,lsp),ppsi(3,lsp),rpsi(3,lsp),torpsi(3,lsp),&
bsp(9,lsp,lst),xsp(9,lsp,lst),zsp(9,lsp,lst),gsp(9,lsp,lst),&
rd(9,lsp,lst),nu(9,lsp,lst),dl(9,lsp,lst),spcos(3,lst),spsin(3,lst),&
rgpsi(3,lsp),cpsi(3,lsp),psirg(3,lsp),psitor(3,lsp),nfpp(3,lsp),&
tfpp(3,lsp))
if(mype==0)then
if(numereq==0)then ! use EFIT & TRANSP data
! EFIT spdata used in GTC: 2D spline b,x,z,rd; 1D spline qpsi,gpsi,torpsi
call spdata(iformat,isp)
! TRANSP prodata used in GTC: tipp,tepp,nepp,zeff,ropp,erpp
call prodata
else
! Construct analytic equilibrium and profile assuming high aspect-ratio,
concentric cross-section
call analyticeq
endif
March 5-7, 2011
GTC Training Course - PKU
26
setup.F90
• Read input parameters
• Initialize quantities determined on equilibrium mesh
!$omp parallel do private(i,pdum,isp,dpx,dp2)
do i=0,mpsi
pdum=psimesh(i)
isp=max(1,min(lsp-1,ceiling(pdum*spdpsi_inv)))
dpx=pdum-spdpsi*real(isp-1)
dp2=dpx*dpx
! 1D spline in psi
meshni(i)=nipp(1,isp)+nipp(2,isp)*dpx+nipp(3,isp)*dp2
kapani(i)=
0.0-(nipp(2,isp)+2.0*nipp(3,isp)*dpx)/meshni(i)
...
enddo
• Allocate fields
! allocate memory
allocate(phi(0:1,mgrid), gradphi(3,0:1,mgrid),...,STAT=mtest)
March 5-7, 2011
GTC Training Course - PKU
27
loadi.F90
• Loads ion markers
!iload=0: ideal MHD with nonuniform pressure profile
!iload=1: Marker uniform temperature & density: n0,te0,ti0 (two-scale
expansion)
!iload=2: Marker non-uniform temperature & uniform density:
n0,te(r),ti(r), marker weight corrected
!iload=3: Physical PDF for marker with non-uniform temperature & density:
n(r),te(r),ti(r)
!iload>99: non-perturbative (full-f) simulation
• Marker data array:
zion(nparam,mimax)
nparam=6 (8 if particle tracking is on)
zion(1,m) –
zion(2,m) –
zion(3,m) –
zion(4,m) – || (parallel velocity)
zion(5,m) – wi (particle weight)
zion(6,m) – (magnetic moment)
March 5-7, 2011
GTC Training Course - PKU
28