300 likes | 549 Views
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
E N D
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 GTC Training Course - PKU
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) GTC Training Course - PKU
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 • Additionally • Magnetic field inhomogeneity is small: GTC Training Course - PKU
Charged particle motion in constant uniform magnetic filed GTC Training Course - PKU
Guiding center motion in the inhomogeneous magnetic filed • Magnetic drift velocity • Magnetic curvature drift • Magnetic gradient drift GTC Training Course - PKU
Gyrocenter equations of motion • Gyrocenter equations of motion [Brizard and Hahm, RMP‘07]: • Modified magnetic field: • Magnetic field perturbation: • ExB drift velocity: GTC Training Course - PKU
Gyrokinetic Vlasov-Maxwell’s system of equations • Gyrokinetic Vlasov equation • Gyrokinetic Poisson’s equation • Gyrokinetic Ampere’s law • Charge density and current GTC Training Course - PKU
Monte-Carlo sampling of phase space volume Fluid moments and potentialsare calculated on stationary grid 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 GTC Training Course - PKU
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 GTC Training Course - PKU
Fluid-kinetic hybrid electron model - I • Simulations of various Alfven modes • Electron continuity equation • Perturbed diamagnetic drift velocity • Perturbed pressure GTC Training Course - PKU
Fluid-kinetic hybrid electron model - II • Lowest order solution of electron DKE based on w/k||v|| expansion • Electron adiabatic response • Inverse gyrokinetic Ampere’s law • Faraday’s law (linear tearing mode is removed) GTC Training Course - PKU
Fluid-kinetic hybrid electron model - III • High-order electron kinetic response • Zonal flow – nonlinearly self-generated large scale ExB flow with electrostatic potential having k||=0 GTC Training Course - PKU
Magnetic flux surface Tokamak geometry • Equilibrium magnetic field (covariant representation) in flux coordinates • Helicity of magnetic fieldis determined by safety factor q (number of toroidal turns of magnetic field line per one poloidal turn) GTC Training Course - PKU
Magnetic coordinates - definition • Covariant representation of the magnetic field • Contravariant representation of the magnetic field • Jacobian GTC Training Course - PKU
Magnetic coordinates - example • Parallel derivative GTC Training Course - PKU
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 GTC Training Course - PKU
Simulation grid • Radial grid in y (i=0:mpsi) • Poloidal angle (q) 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 GTC Training Course - PKU
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 Core 0 Core 1 Core 0 Core 1 Core 2 Core 3 Core 2 Core 3 Node 0 Node 1 GTC Training Course - PKU
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,operation,comm,ierror) GTC Training Course - PKU
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 GTC Training Course - PKU
GTC structure dne dge1&dfi dA|| Dynamics due dA|| ZF dfes dfind Fields dA|| dui dne1 dni due1 dne Sources GTC Training Course - PKU
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 GTC Training Course - PKU
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) GTC Training Course - PKU
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 GTC Training Course - PKU
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 GTC Training Course - PKU
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) GTC Training Course - PKU
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) – y zion(2,m) – q zion(3,m) – z zion(4,m) – r|| (parallel velocity) zion(5,m) – wi (particle weight) zion(6,m) – m (magnetic moment) GTC Training Course - PKU