Home

HMGC User guide

1. uu 0 f0 gt Geu t _f0 x2 397E 03 Ea O fO x2 911 053 z a 1 9 1 oo 3 3 g 3 3 1 sn Be 2 n SM E g o5 LL o5 LL o5 LL ga sa sa ve di 2 ve di 2 ve di Q 9 T 9 9 a2 caz a2 uy 1 S 1 B Es ut ug a a a 24 a amp O cas O cas O E ET ET 4 0 2 0 0 5 1 0 5 1 D 025 Db 0 75 1 u l E E Figure 14 Left fg r u u Centre fy r u u transformed on the E a mesh Right fy r E a All figures are obtained summing from rjz9 p mi 0 375 to Tjero pi max 0 625 The relaxed time is frelax 120 With respect to the in put file xplot_deltaf_input shown in the text the figures are obtained by assigning iflag df 1 All figures refer to the frame 211 two 126 52 5 5 Output file power data The file power data contains the time series of the power exchanged between the energetic particles and the waves Py r u thus integrated over 0 p and is read by the program plot power f The power is stored on a r u u grid defined for plotting purposes The file is in unformatted format The data stored on the 4 u mesh can be transformed and plotted on a E a one note that this is not equivalent to store directly the distribution function in the E a space indeed the single particle energy dependence on is lost The data can be also summed up between to radii to reduce particle noise see Fig 15 An example of the input file xplot_power_input f
2. write 64 Capmhdg im jer psieqg im jer im 1 1mpert amp jer 0 nlr Note that the file APWRITE contains the Ym n r components perturbation part apmhdg plus equilibrium part psieqg and the file PHIWRITE contains 5 r perturbation part phmhdg equilibrium part assumed to be zero ERG is radial coordinate of the mesh used by the the gyrokinetic module LMPERT LM Those files are used to produce a sequence of frames of a series of quantities only the ones produced for m n r are shown in the following Figs 8 10 using the plotting program plot field f The plots are the following con tour plot of r 0 phiplot r 0 yi t with superimposed the trajectory of the i test particle trajectories of the i test particle r r t 0 0 1 pi yilt uj u t see Fig 8 trajectory of the i test particle in the poloidal R Z and equatorial X Y plane see Fig 9 radial profiles of the m n Fourier components contour plot of m n r in the plane r m with superimposed the curve m nq r contour plot of the power spectrum P r w in the plane r w with superimposed the lower and upper Alfv n continua for the toroidal gap r radial component of the displacement and Te r electron temperature fluctuation assuming incompressible perturbations The power spectrum P r w is defined as P r w UmnPmm 74 X Umyn Im n r w T l m n r w 10 45 An example of
3. 0 call cerca_massimi for plotta_max 0 no 1 yes 0 iflag_deltate synthetic diagnostic Delta_Te 0 1 0 1 csi0 csil min max in csi if 0 O use max available interval 50 50 0 deltate0 deltate1 min max in deltate if 0 O use max available interval Note that the previous input file will produce a sequence of plots starting from twao 0 to twao 96 only a single frame will be shown in the following Figures 8 9 10 iflag test can vary from 0 no test particle plots to NPTEST producing plots for the i th test particle iflag trajectory 1 will produce in addition the trajectories of the selected test particle rtest t Atest t Ptest t Utest t The input file Te vs erg DIII D 1 interp txt contains the vectors r T r on the NRL mesh for synthetic diagnostics purposes to be produced by the user Note that for producing the power spectrum in the plane r w one has to choose the time window used in the FFT twindow Tpprwyo the minimum frequency wmin resolved is 46 given by wmin Wao 27 twindow whereas the maximum frequency is given by Wmax Wao Tt At v NWRITE DT To minimize the effect of having a finite time sequence the data can be multiplied by a Hanning function which essentially is a function picked at the middle of the time sequence and which goes to zero toward the edge To increase the number of points in the frequency direction but not the content of information a buffer
4. 8 x 21 168 has been used The factor 2 in the variable nintphi nintphi 2 nmax 1 and the value of 4 for the variable nph_su_nintphi are such that NPH is NPH 8nmax 8 2 16 Those values are the ones typically used in HMGC simulations The quantity nne should be such that npart nne nnalpha nlr nth nph ippc 3 a simple program to compute the optimal values to distribute the particles in the E a space is provided calcolo_nne f for given nlr nth maximum n nph_su_nintphi ippc The program asks at the beginning which source you are referring to enter 0 for data referring to e3_complete 24 Quantities Comments NLR number of radial cells for the gyrokinetic module NLR 1 mesh points NTH number of points in the poloidal angle direction nintphi nintphi 2 nmax 1 Parameter for y toroidal angle mesh the energetic particle distribution function will be loaded on nph_su nintphi toroidal mesh points and then replicated on nintphi 1 remaining sectors nph_su_nintphi parameter for y toroidal angle mesh NPH NPH nph_su_nintphi nintphi number of points in the y toroidal angle direction ippc number of particles per cell in each direction of the physical space r 9 p nne number of particles in the energy space nne nnalpha npart nintphi nnalpha is the number of particles in the pitch angle direction npart nlr nth nph ippc 3 total number of particles
5. Ala B pa L pa 1 r9 pa BETA RHOA AGROWTH parameters for an ad hoc growing term added BGROWTH in the vorticity equation to simulate some parti CGROWTH cle drive for CGROWTH lt T lt CGROWTH DGROWTH DGROWTH GROWTHR I GROWTH EXP RCI AGROWTH 2 BGROWTH else GROWTHR I 0 now obsolete ITAERSP ITAERSP 1 should drive a 1 1 2 1 TAE using the ad hoc driving term ITAERSP 2 should drive a 1 1 2 1 RSP Resistive Shear Periodic mode requires GROWTH ne 0 0 and n 1 m 1 2 modes SMOFAC amplitude of a smoothing factor to control numerical instabilities in the center r 0 hyper resistivity and hyper viscosity terms AMP complex amplitude factor for the initial perturbation GROWTH amplitude of the ad hoc growing factor GYRO logical variable ifGYRO true the gyrokinetic module is called CYLIN logical variable ifCYLIN true MHD module considers cylindrical limit while gyrokinetic module retains finite correction BISEC logical variable ifBISEC true allow time bisection in the MHD mod ule and hence in the gyrokinetic one SKIN electron skin depth if 0 electron inertia neglected now obsolete not used epsill parameter used to reduce e correction at the edge occasionally used to control edge numerical instabilities arising from MHD module see Eq 8 cgeps parameter used to reduce e correction at the edge see Eq 8 cweps parameter us
6. Tar Electron temperature eV Bu r Beam beta toroidal dimensionless Note that this transp value is reduced by fast ion transport wr Toroidal angular velocity rad s q r safety factor Table 10 Experimental radial profiles provided by DIII D team The quantity By r wil not be used as an input data but only to compare with the computed y r from the code this latter quantity will depend indeed on the model used to load the initial distribution function The quantity w r could be used to compare the experimental frequency vexp spectra with the ones obtained by the code Mode vexp Vcode n 27 wy with n the toroidal mode number In Fig 18 the geometry of the neutral beam is shown The beam is essentially tangential wit injection angle ag arccos Rian Ro arccos 1 15m 1 688m 47 055 deg The average birth energy is Eg 0 077 MeV 2 3 x 0 075 MeV 1 3 x 0 081 MeV From the data profiles provided and from global information a single toroidal mode number n 2 simulation can be set up performing the following steps 58 SI i Electron temperature eV Lon temperature eV Impurity density carbon cm 3 B Beam ion density cm 3 Note that this transp value is 3 310 7 7 10 vlad GV 1 vlad GV us ape d 110 i 210 J foe Po o 3 e m dereen suo g 2108 TS s 3 oe 5 L da m g 810 F o J
7. NMODOM number of Fourier components for the gyrokinetic module they must be the Fourier modes of the MHD module plus two poloidal satellite modes for each toroidal mode considered in the simulation NRZ the particles are evolved in a NRZ NRZ grid in the R Z plane around the r 0 point to avoid problems related to the singular point r 0 in polar coordinates Table 5 Structure of the file grid inc 25 3 7 Input file KININP This file is the main parameter input for the gyrokinetic module Here is a sample for the DIII D simulation see also Tables 6 7 0 032863457d0 RHOSA idistr 1 sqrt T_HO m_H Omega_cH0 a idistr 2 sqrt E_0 m_H Omega_cH0 a idistr 3 sqrt T_perp_H0 m_H 0mega_cH0 a 0 271063836d0 VTHSVA idistr 1 sqrt T H0 m H v AO idistr 2 sqrt E 0 m H v AO idistr 3 sqrt T perp HO m H v AO 1 5d0 sigma O0 idistr 3 T perp HO T par HO 2 3256d0 usdelta_input idistr 1 2 anisotropy parameter 1 width 0 68128d0 cosalfa_0_input idistr 1 2 anisotropy parameter cosine of injection pitch angle 4 153850158d0 e0sec0 idistr 2 on axis EO E critO 2D00 ALF parameter controlling non uniform radial particle loading 0 NPIC parameter controlling non uniform radial particle loading if npic ne 0 erO i del i i 1 npic must follow 0 264848976d0 ENHSNI n_HO n_iO on axis ratio between energetic particle and bulk ion densities 1 D00 EMHSMI m H m i ratio between energetic particle and bulk
8. lar component of the wave vector to the magnetic field and py the energetic ion Larmor radius though fully retaining magnetic drift orbit widths and solved by particle in cell PIC techniques The coupling between energetic ions and thermal plasma is obtained through the divergence of the energetic ion pressure tensor which enters the vorticity equation Numerical simulations of experimental conditions are performed by fitting the relevant thermal plasma quantities the on axis equilibrium magnetic field major and minor radii Ro and a respectively the safety factor q the electron n and ion n plasma densities the electron temperature T the anisotropy of the energetic ion distribution function and the ratio Gy between fast ion and magnetic pressures In order to retain the relevant finite Larmor radius effects without resolving the details of the gyromotion the energetic ions are evolved in their gyrocenter coordinate system which corresponds to averaging the single particle equations of motion over the fast Larmor precession HMGC has been successfully applied to the interpretation of the experimental evi dences of rapid transport of energetic ions related with fluctuations in the Alfv n mode frequency range in auxiliary heated JT 60U discharges in connection with so called Abrupt Large amplitude Events ALEs 3 4 5 HMGC results have also suggested a possible justification of the large discrepancy observed in reversed sh
9. The output of the program eqe3aaab is a file named EQNEW This file will be copied by the HMGC execution script to the file named INCOND Its structure is shown in table 2 12 Quantities Comments the EQUIPA namelist see Table 1 0 DO NR 1 0 FORMAT 3120 NR is the number of radial grid points R I I 1 NR FORMAT 2D30 15 the normalized radial coordi nate a sequence of radial profiles for the m 0 n 0 and m 1 n 0 Fourier components for w and resistivity profile 7 in the following format two blank lines PSI 1 a line containing the following text PSI 1 fotu alt PHI 3 for m n r or RES 4 for Nmn r real m real n FORMAT 2F20 0 m n being the poloidal and toroidal mode numbers respectively for the equi librium is n 0 real Wmn I imag m n I FORMAT 2E30 15 I 1 NR 1 only NR points for and n Table 2 Structure of the file EQNEW 13 Note that here n 0 equilibrium also note that the Ym n r harmonics have one more radial point NR 1 corresponding to the position of a resistive wall this option is usually not considered The electrostatic scalar potential dm n r components for the equilibrium are identically zero equilibrium without fluid flow and usually MHD module of HMGC used in linear mode the resistivity profile is taken constant in radius 70 0 1 71 0 0 Note also that HMGC defines the vo to be Yo r
10. Z 0 top view of the torus dotted line refers to r 0 the cross indicates the axis of symmetry of the torus produced by assigning iflag test 1 and iflag trajectory 1 P Hio 3x 576 0 tow 288 00 Ema A IN rs m ley n 2 tww 96 00 aupra wr sia Moda 9 255 416F 10 4 847E 06 110 3 2 9 o 2 52 5 a a 0 8 Pic 9 al 9 2 2 e x 02 c m 0 6 15 2 lt lt 17 2 T 19 2L LJ J LJ 0 4 ete Q 14 o o o ua a ua 0 2 8 8 i 8 Q Q 1 Q C H 9 4 4 5 0 z 0 0 20 2 04 06 08 1 0 0 2 04 0 6 08 1 0 20 2 04 0 6 08 1 ria ria ria Figure 10 Left mn r produced by assigning iflag fourier comp 1 and iflag_contour 0 Centre mn r contour plot and superimposed the curve m nq r Right frequency spectra of mn in the plane r w contour plot with superim posed the upper and lower Alfv n continua of the toroidal gap for a particular toroidal mode number automatically chosen from the first perturbed mode in this case n 2 produced by assigning itot 961 ipl0 ipli 481 iflag power spectrum 1 and other data from file xplot_field_input 48 f Gap t n Ato 1x 48 0 tou 144 00 gt o c th ar Ato 1x 48 0 to 144 00 gt ero fh arr Ato 3x 48 0 toe 144 00 Sum over m n I min l max 3 23 Sum over m n I min l abet 125 Sum over m n l min l max 3 23 x 2 14 Li 0026 13 x DOBE 13 fi fi 1 BE 12 HMGC code ENEA Frascati HMGC
11. ibuffer 3 i e din Onan Oe erus e Gets ees up UM SIEHE S Frequency spectra during the saturated phase twao 432 and different values for twindow first row twindow 72 Wmin Wao 0 0873 second row twindow 144 Wmin Wao 0 0436 third row twindow 288 Wmin wao 0 0218 ihann and ibuffer first column ihann 0 ibuffer 1 second column ihann 1 ibuffer 1 third col umn ihann 1 ibuffer 3 i e Wmin plot Wmin 9s adi oe Left fy r p u Centre fy r u u transformed on the E a mesh Right fg r E o All figures are obtained summing from TjerOplmin 0 375 tO Tjer0plmar 0 625 Note that the input file xplot_deltaf_input shown in the text refers to the left plot the fig ure in the centre is obtained with iflag_versus 2 the figure on the right is obtained reading the file deltaf ealfa data All figures refer to th frame Paluan s G Lacer dto eA a HEU ee di D AASA oe ndr 48 14 15 16 17 18 Left dfy r u u Centre 6fy r u u transformed on the E a mesh Right 6fy r E a All figures are obtained summing from rjero_pi_min 0 375 to rjero_pimar 0 625 The relaxed time is treiax 120 With respect to the input file xplot_deltaf_input shown in the text the fig ures are obtained by assigning iflag df 1 All figures refer to the frame A TIZIA RENO ETERO Left Pg r u u Also the curves representing the maximum energy loaded in the initial distribution function dotted lines and the trapped unt
12. the n 2 reference DIII D case right a n 2 4 case Black dots represent the modes actually included in the simulations the red crosses represent the modes considered in the simulation because of complex conjugate condition 16 3 2 Input file PARAM This input file contains the main input parameters for the MHD module and for some general input parameters for the simulation Here is a sample for the DIII D simulation see also Tables 3 4 160 3 20 92 1 0D 5 1 0D 8 0 02 30 1 000001 1 0D10 0 0D00 0 0D00 TRUE TRUE FALSE FALSE 9173d0 69776d0 6471d0 T 01 05 95 FPRrROODOODOOCOOWF DO 1 D 07 1 0 00D 0 TRUE FALSE FALSE 0 D 2 0 10d0 0 95d0 0 025d0 DIAPOS amp END NCYCLE Number of GK calls Number of Time Steps NTS NCYCLE NSUBCY NSUBCY Number of MHD calls per each GK call NOUT Number of time sequences total time steps NTS NOUT LR1 Maximum value should be LR1 4 LM ETA standard value 1 0d 5 VISC standard value 1 0d 8 DT standard value 0 02 NPRI NTS NPRI lt MAXPRI RWALL resistive wall normalized radius parameter required but not used by fortran give any real number TAUWAL resistive wall characteristic time parameter required but not used by fortran give any real number VEDGE plasma bulk velocity at the edge parameter required but not used by fortran give any real number CURAMP current ramp ios ignored FREZOO 1 1 m n
13. 1 iprek 1 iprek 1 0 1 th Fourier component of the PREK off iprek 1 1 on 21 1 iprek 1 26 Quantities Comments RHOSA pu a particle Larmor radius normalized to minor radius com puted 1 at the on axis energetic particle temperature Tyo py a VTuo ma Qero a if a Maxwellian distribution is considered IDISTR 1 2 at the birth energy Eo px a VEo my Qe40 a if a slowing down distribution function is assumed IDISTR 2 3 at the on axis perpendicular energetic particle temperature Tino pufa VT u0 mg Qen0 a if a bi Maxwellian distri bution is considered IDISTR 3 VTHSVA Uin VaAo ratio between energetic particle thermal velocity and on axis Alfv n velocity for IDISTR 1 vn v4o V THo Mu Va0 for IDISTR 2 v n v40 V Eo mu Va0 for IDISTR 3 vn v40 VTL Ho MH VA0 sigma_0 ratio between on axis perpendicular and parallel energetic particle temperatures used only if IDISTR 3 usdelta_input parameter for anisotropic particle distribution function it corre sponds to the inverse of the width A see Sects 6 1 6 2 of the distribution function around the injection pitch angle ap for Maxwellian or slowing down distribution functions IDISTR 1 or IDISTR 2 cosalfa_0_input COS Qo cosine of the injection pitch angle ag for Maxwellian or slow ing down distribution functions IDISTR 1 or IDISTR 2 e0sec0 Eo Ecrito on axis ratio between birth energy and crit
14. can be derived from the radial profiles of the above listed plasma quantities e g using a spreadsheet and it will be used as temp exp DIII D 1 file table of v r Eci r Ecrit 61 Similarly the energetic particle normalized density profile will be derived by the Beam ion density profile and used as den_exp_DIII_D_1 file table of v r ng r ngo For each of the previous files the program interp spline f is called by the execution script producing the files den spli data and temp spli data 7 5 How to setup an HMGC run energetic particle dimensioning file grid inc The parameters which describe the energetic particle dimensioning have been chosen according to Table 11 NLR 64 low n should be such to resolve the typical radial width structure of the perturbed mode NTH 168 usually take 8 poloidal mesh points times Mmax nintphi 2 x n 1 nph su nintphi 4 ippc 2 two particle per cell per each direction r 0 nne 672 as obtained by the program calcolo nne f NMODOM 2T see Table 5 NRZ 5 it depends from the Shafranov shift of the equi librium see output of the program eqe3aaab f typically NRZ such to include the magnetic axis is taken Table 11 Preparing the grid inc file 7 6 How to setup an HMGC run energetic particle parameters file KININP The energetic particle parameters can be calculated from the experimental data provided namely energetic io
15. code ENEA Frascati HMGC code ENEA Frascati 0 2 0 4 06 08 1 0 02 04 06 08 1 02 04 06 08 1 r a r a r a f Gap keto rP Bto 1x 96 0 tww 144 00 gt e g Iplar Atos 1x 98 0 to 144 00 gt ceo i ni Ato 3x 96 0 to 144 00 Sum over m n I T Ippica 3 23 Sum over m n I min l max 3 23 Sum over m n l min l max 3 23 X SE TO aa X SAT ESTA x 6315F TS qaaa i 1 dad 0 44 0 24 0 24 HMGC code ENEA Frascati HMGC code ENEA Frascati HMGC code ENEA Frascati 0 02 04 06 08 1 0 02 04 06 0 8 1 0 02 04 06 08 1 r a r a r a W Gap Ipla rj Ato 1x 144 0 to 144 00 f Gap Iplar Atos 1x 144 0 tow d hii 00 G9 an Iplar Ato a 3x 144 0 to 144 00 Sum over m n I min z max Sum over m n l min RIE 3 23 Sum over m n I min l mas 25 _ _ PE are x D 2E US i PARE PI 935E7 i i A prodi BE 12 E Q Q Q A a A a a g g g LL LL LL 0 44 21 Mp 27 x x x od Lu LJ LJ 0 25 o o 7J 3 3 a 9 a o o o e e 9 0 2 04 06 08 1 F 0 02 04 06 08 1 F 0 02 04 06 08 1 F r a ria r a wf Gap Iplar Aton 1x 192 0 top 144 00 o f eg ar Ato 1x 192 0 to 144 00 f Gap Iplar Aton 3x 192 0 teyo 144 00 Sum over m n l min D TIE US Sum over m n I min l reba 3 23 Sum over m n l min ose 3 23 FER IE Lu 1 fior x 2 4 28E U8 1 1 Lu 92F 10 1 Q Q Q 9 A a A J 9 A 2 2 2 0 44 v o
16. down definition of the distribution function it follows that sinada0 a 2 19 0 54 Relevant limits are A O a 1 20 A 0 O a 26 cos a cos ag 21 The definition of the pressure tensor is pj fato 2FEcosa 22 p f d v F Esin a 23 2 p iugis 3 Pure 24 where F is the distribution function From the normalization condition 19 it follows PI PN 25 pi SpA 26 ire tas Bee 27 hen sf dus aula 28 AL A Ss db 29 Then pi pj E 30 If the ratio p pj is provided experimentally the corresponding value for A can be obtained implicitly from Eq 30 see Fig 16 The explicit expression for Aj is Ay f dasina cos aola 31 2A24 1 2 erf i cosan erf teosan x Cm ray CER 95 2 cos ag RA Pperpsuppar Delta o a cosalfa_0 0 68128 Figure 16 Pperp Ppar VS A 6 2 Slowing down distribution function idistr 2 The slowing down energetic particle equilibrium distribution function fgo is defined as follows n 3 MH 3 2 JEQsa TA gt JEQ sd 32 Ecrit 0 fi O a ap A Han 4 a ao A 33 TUP E Eere 24 1 In E Eq Ernie J In this case the quantity 7 is given by 7 Ecrit 4 Ecrito with Ecrit being the critical energy see Stix 8 2 3 A nlp Z2 Es 6 MeV 14 87 v MeV i 34 with Ay the mass number of the energetic ions A and Z the mass number and ele
17. ion masses 999 99d9 timkin_anu energetic particle density ramps for t gt timkin_anu set to a large value to avoid ramping 120 d0 TIMKIN_RELAX time at which the distribution function is assumed to be relaxed used for ramping and diagnostics 1 166212d0 ANU_MAX energetic particle density ramping parameter 0 166212d 2 ANU_DOT_O energetic particle density ramping parameter 0 i_write_deltaf 0 no output 1 output of f r mu u t 2 output of f r E alpha t 0 i write power 0 no output 1 output of wave particle power exchange P r mu u t 30 NWRITE energetic particle quantities written every NWRITE time steps 2 IDISTR 1 Maxwellian 2 Slowing down 3 bi Maxwellian 1 IDELTF O full f 1 delta f 0 ILIN 0 fully non linear GK simulation 1 linear GK simulation 0 IRLSRO 0 R_1 with epsilon corrections 1 R l1 RO 1 IMIRR 0 mirroring term off 1 mirroring term on 1 IWOO 0 grad B drift contribution to the source term neglected 1 retained 1 ILANDA 0 Landau damping term off i Landau damping term on 1 ICURV 0 curvature term off i curvature term on 1 IOMST 0 omega star term off 1 omega star term on 2 NPTEST number of test particles 0 ITEST 0 true test particle 1 simulation particle 5 0 0 1 1 EROT THOT PHOT AMOT UOT for t 0 r theta phi mu u 1 ITEST 0 true test particle 1 simulation particle 1352 LTEST simulation particle identification number 0 ER_PERTO energetic particle pressure term PREK set to zero for r lt er_pert0 11
18. is the equation of state and generally describes the polytropic evolution of the plasma It may be combined with the continuity equation and written as TPV v 47 Equation 42 the so called ideal Ohm s law describes the plasma as a perfectly conducting fluid from which the expression ideal MHD originates Note that in the more general case 66 in which plasma resistivity n is considered eq 42 is replaced by 1 E vxB nJ 48 The elliptic operator A is given by Vw O 10v Oy A 4 RV 4 tay Hop 755 an i 9 2 Order reduced MHD We introduce here a simplified version of the resistive MHD equations which will be useful later for the non linear study of the Alfv n modes greatly reducing the complexity of the problem We start from the resistive MHD equations 39 45 with the Ohm s law generalized by eq 48 Since tokamak plasmas are characterized by values of the safety factor g r rBy RBs O 1 By and By are respectively the toroidal and poloidal component of the magnetic field and inverse aspect ratio e a Ro much lower than unity MHD equations can be simplified by expanding in powers of e This procedure has been widely used since the first paper of Strauss 11 both for analytical and numerical work At the leading order in e O c and considering the low approximation 8 O c the reduced MHD equations describe the plasma in the cylindrical approximation
19. is useful in reducing numerical noise as long as F y F pol In the present paper when adopting the 6F approach we take F yo to be Maxwellian QyM e DI F uo Ce ng R exp Ty where ny R and Ty are respectively the energetic particle equilibrium density and uni form temperature The r h s of eq 77 is then given by S t R M U Fro Pci i vie dt Ty mg On MH egU a a CHE epu RUE Vox Vaj 82 We also assume the following model for the energetic particle equilibrium density nur nan exp 5 83 74 where nyo is the on axis density End of excerpt form the Appendices of Ref 9 References 1 Frieman E and Chen L 1982 Physics of Fluids 25 502 2 Lee W 1987 Journal of Computational Physics 72 243 3 Shinohara K et al 2001 Nuclear Fusion 41 603 4 Shinohara K et al 2002 Nuclear Fusion 42 942 5 Shinohara K et al 2004 Plasma Phys Contr Fusion 46 S31 45 6 G Vlad S Briguglio G Fogaccia and B Di Martino Gridless finite size particle plasma simulation Computer Physics Communication 134 58 77 2001 7 G Vlad S Briguglio G Fogaccia F Zonca W W Heidbrink and M A Van Zeeland Particle simulations of Alfv n modes in reversed shear DIII D discharges heated by neutral beams In 10th IAEA Technical Meeting on Energetic Particles in Magnetic Confinement Systems Kloster Seeon 8 10 Oct 2007 page OT 7 Vienna Austria 2007 Intern
20. of motion eqs 65 Such markers can then be interpreted as the phase space coordinates of a set of Npart particles and G t Z can be approximated by Npart G t Z Y A t G t Zi t 8 Z Zi 68 l 1 12 The time variation of the volume element A t is then given by dA 8 dZ e At Z 69 OZ t Z t For the purpose of calculating the pressure tensor components eq 63 it is convenient to directly represent the quantity D zF g according to its discretized form Npart D z t Z Fu t Z M w t 6 Z Z 70 l 1 with the weight factor w defined by and A Ai t D z t Zi t 72 In fact from eqs 64 69 and from the Liouville theorem a dZ D l D z2 t 2 gt Z az zenZ dt 0 73 it is immediate to show that dA an 74 and a D At each time step the fluctuating electromagnetic potentials are computed at the grid points of a three dimensional toroidal domain in terms of the Fourier components yielded by the field solver Phase space coordinates are then evolved in the fluctuating fields and the pressure tensor components at the grid points are updated in order to close the MHD equations for the next time step Field values at each particle position are obtained by trilinear interpolation of the fields at the vertices of the cell the particle belongs to The corresponding trilinear weight function is adopted after pushing the parti
21. produces the output of the distribution function i_write_deltaf 0 no output is produced iwrite_deltaf 1 f r u u t i write deltaf 2 f r E a t u u E and a are the magnetic moment parallel velocity particle energy and pitch angle respectively i write power flag for writing the power exchange P r u u t between particles and waves 0 no output is produced 1 output is produced NWRITE energetic particle quantities are written on output files every NWRITE NSUBCY time steps IDISTR 1 Maxwellian 2 Slowing down 3 bi Maxwellian IDELTF 0 performs a full f simulation 1 performs a f simulation standard use of HMGC Table 7 Structure of the file KININP continued 28 Quantities Comments ILIN 0 performs a gyrokinetic non linear simulation standard use of HMGC 1 performs a gyrokinetic linear simulation note that ILIN 1 corresponds exactly to a linear simulation only if the initial distribution function is a true equilibrium function IRLSRO for the generic simulation particle 0 retains e corrections R Ro 1 e cos 0j standard use of HMGC 1 approximates R Ro IMIRR 0 mirroring terms off 1 mirroring terms on standard use of HMGC IWOO IW00 0 causes to neglect the contribution of the grad B drift to the source term in the delta f Vlasov equation 1 terms are retained standard use of HMGC ILANDA 0 Landau damping off 1 Lan
22. 0 0 FREZ10 1 2 m n 1 0 EQUIL DROP NPROFI 0 DEN 1 1 DEN DEN RHOA ALFA BETA 2 DEN Q Q0 2 ALFA BETA RHOA AGROWTH ad hoc growing factor parameter BGROWTH ad hoc growing factor parameter CGROWTH ad hoc growing factor parameter DGROWTH ad hoc growing factor parameter ITAERSP 1 drives TAE 2 drives RSP requires GROWTH ne 0 SMOFAC amplitude of the smoothing factor D 07 AMP complex amplitude factor for the initial perturbation GROWTH amplitude of the ad hoc growing term GYRO call gyrokinetic module CYLIN true MHD cylindrical limit BISEC true bisection allowed SKIN el skin depth skin 0 DO gt el inertia neglected epsill parameter used to reduce toroidal corrections at the edge cgeps parameter used to reduce toroidal corrections at the edge cweps parameter used to reduce toroidal corrections at the edge NRCHNL 6 RCHNL 1 0 200 RCHNL 2 0 300 RCHNL 3 0 400 RCHNL 4 0 500 RCHNL 5 0 650 RCHNL 6 0 800 17 Quantities Comments NCYCLE number of calls of the gyrokinetic module for each of the NOUT time sequences NSUBCY Number of MHD calls per each gyrokinetic call NCYCLE NSUBCY is the number of calls of the MHD module for each of the NOUT time sequences NOUT number of time sequences total time steps NCYCLE NSUBCY NOUT LR1 maximum number of modes per MHD fields which are read in the file INCOND maximum value should be LR1 4 LM ETA normal
23. 1 1 DT NWRITE 96 The parameter 5 in the last line of the file xplot density input is a parameter required by the plotting routines HIGZ from CERN which identifies the graphic window The plotted quantities are given in terms of the quantities written in the file OUTDNT by e npg r ngo DNTOT r 27Ar where Ar 1 NLR e Gi r 2 x EMHSMI x VTHSVA x ENHSNI x 2 x PPERP r 3 PPARA r 3 27Ar e ay Ro a x q r dBy dr Al Nu niu tw 96 00 Bur tw 96 00 eur twy 96 00 zia bacca tipi di 1 tii zia db oaa ili a db aca bau BL lr aaa aca 0 01 HMGC code ENEA Frascati HMGC code ENEA Frascati HMGC code ENEA Frascati 0 42 0 4 0 6 08 1 0 0 2 0 4 0 6 0 8 1 0 02 04 0 6 0 8 1 r a ria ria Figure 7 Example for density 9g red dashed curve is 3 7 green dotted curve is 8 and ay radial profiles 42 5 2 Output file TESTWRIK This file contains the time history of the test particles and it is used to plot test particle quantities see Fig 8 using the plot program plot field f The quantities written are Trest test test Utest Wtest Atest see Sect 9 4 The following excerpt shows the write state ments in the file TESTWRIK 1 IF ISTEPO EQ 0 THEN WRITE 29 ASPECT WRITE 29 NPTEST DO 335 I 1 NPTEST WRITE 29 ITEST I IF ITEST I EQ O THEN WRITE 29 EROT I THOT I PHOT I AMOT I UOT I ELSE WRITE
24. 1 0 and positive in the plasma 14 3 Execution script The execution script of HMGC prepares a number of files used for compiling and running HMCG Hereafter is a list of them referring to a DIII D case see fig 2 simulation which consider an equilibrium with m 0 1 and n 0 modes this is mandatory and a perturbed n 2 mode with poloidal components ranging from m 1 to m 21 Note that because of symmetry conditions in the Fourier space only modes in a half m n plane are required the other ones being considered using the rule Wmn T v r reality of 7 0 The choice of considering only the mode in the positive half plane defined by n gt 0 is used More over the conventions for the Fourier transform are v r 0 p Vo o r zs 6 2 p Re Um n r cos m0 ng IM Umn r sin m ny 1 2 LM Uma apy LD MTA 7 j 1 No k 1 No with l being the mode index m m l n n l LM the total number of Fourier com ponents in the simulation see Sect 3 4 Ng and N the mesh points of the 0 and g grids respectively The choice for the poloidal Fourier components included in the simulation derives usually from considering Mmin Nqmin Mmax NGmax Some restrictions could be imposed by FFT requirements see Sect 3 1 15 3 1 Include file modi inc Parameter definitions for compiling the MHD module e3 complete F of HMGC NR is the MHD radial grid must be NR NREQ see Table
25. 1 LM is the number of Fourier components considered in the simulations NMAX Nmax 1 is the maximum toroidal mode included in the simulation n plus unity MMAX is the maximum number of poloidal Fourier components for fixed toroidal mode number n MAXPRI is a parameter to dimension a buffer for certain output quantities In Fig 3 are shown two sketches of the m n plane used by HMGC for better clarify the parameters meaning A constraint given by FFT routines impose that 4 MMAX 1 is a valid number for the FFT see e g http publib boulder ibm com infocenter clresctr vxrx topic com ibm cluster ess 43 guideref doc am501 formul2 html Actually the IBM ESSL package is used but routines which use NAG modules are also included in the source files although they could be out of date PARAMETER NR 150 NMAX 3 MMAX 21 LM 23 MAXPRI 200 e m_n0 sim e m n2 sim e m n4 sim x m n0 cc X mn o x m n4 cc TMODE 1 TMODE 1 diii ronis TTT ARAARA AAAA MMAX 21 IMMAX IA 1 HH NMAN 241 H 2 BM 799999 z 0 H s 0 1 c c modes not explicitly 2 c c modes not explicitly included iih the simulation included in he simulation 2 4 LIL 5 10 15 20 15 10 5 0 5 10 15 m m Figure 3 Fourier modes included in a HMGC simulation left
26. 1 anisotropic slowing down idistr 2 and bi Maxwellian idistr 3 Their expressions are listed hereafter 6 1 Maxwellian distribution function idistr 1 The Maxwellian energetic particle equilibrium distribution function fEQ Maxw is defined as follows ngo MH 3 2 4 fEQ Maxw TP fEQ Maxw 11 A _ AQ E Tg 4 SEQ Maxw roy 9e ag A e 12 E pun wen 13 COS COS Qr 2 4 exp exeo 60 a ao A 14 6002 AU ef ESS erf Hem us cosa 15 V 2E mg D z pOH sina E 16 nu w nxon y 17 TH w TuH ot 18 with E the energy ng v the radial density profile Tg v the temperature u the parallel to the equilibrium magnetic field velocity u the conserved magnetic moment a the pitch angle of the energetic particles and O a represents the anisotropy of the distribution function ngo ng r 0 Too Ty r 0 Qeg ex B myc with eg my the charge and the mass of the energetic particle respectively and B the local equilibrium magnetic field The parameters entering in the anisotropy term O have the following meaning o is the injection angle e g in the case of beams and A is the width of the beam around cos ao In the code the parallel velocity is normalized to the on axis energetic particle thermal velocity u vmo with viho V Tgo mng and the magnetic moment is normalized to fi uQcHo THo with QeHo eg Bo mye From the slowing
27. 29 LTEST I ENDIF CONTINUE DO 1 L 1 NPTEST WRITE 29 ERTEST L THTEST L PHTEST L UTEST L WRITE 29 WTEST L DTEST L CONTINUE Units and normalizations are ERTEST rtest normalized to minor radius a THTEST brest in radiants PHTEST test in radiants UTEST Utest Upar vn 4 0 with va pb 0 VTa 4 0 mu Tu 4 0 is 1 the temperature at v 0 for Maxwellian distribution function idistr 1 2 the birth energy Eo for the slowing down distribution function idistr 2 3 the perpendicular temperature at i 0 for the bi Maxwellian distribution function idistr 3 AMOT u Og ib 0 TH w 0 with u being the magnetic moment and Ny y 0 the Larmor frequency of the energetic particle at Y 0 43 Note that here p 0 defines the magnetic axis 44 5 3 Outputs file PHIWRITE and APWRITE These files contain the time sequences of the radial profiles of the Fourier components of m n r and Ym n r respectively on the radial grid of the gyrokinetic module and normalized according to the normalizations used therein The files are in unformatted format The following excerpts show the write statements in the file PHIWRITE unit 63 and APWRITE unit 64 WRITE 63 ISTEPO TIMKIN IF ISTEPO EQ 0 THEN WRITE 63 NLR NTH NPH WRITE 63 ERG JER JER O NLR WRITE 63 LMPERT write 63 mmode im nmode im im 1 1mpert ENDIF write 63 phmhdg im jer im 1 l1mpert jer O nlr
28. 3 3 o Bm o E 1108 j 9 610 L o J b 5 E ti amp 410 f E 510 j a re 210 eo d e ki 0 0 1 i i aci 0 0 0 2 0 4 0 6 0 8 1 0 0 0 0 2 0 4 0 6 0 8 1 0 index 19 index 19 Toroidal angular velocity rad s 410 11 T vlad GV ivlad GV 410t J 10 d o E31 E j oL j gt o 310F i 8E E 2 210 o 7 S E E di 1 amp gt 10 F uf 1 6E amp 4 e o 3 110 EF j she i E 3 Tog ao S510 Peo 64 4L tegas J 0 1 1 f f 3 1 1 f 1 0 0 0 2 0 4 0 6 0 8 1 0 0 0 0 2 0 4 0 6 0 8 1 0 index 19 index 19 Figure 17 Experimental profiles for the DIII D shot 122117 59 9 Beam beta toroidal dimensionless Note that this transp value is re vlad GV 0 2 1 1 0 4 0 6 index 19 Rign 115 cm beam Ro 168 8 cm COSO Ryjan Ro tan Figure 18 Beam geometry 7 1 How to setup an HMGC run preparing the equilibrium file EQNEW Using the experimental q r profile a EQNEW file should be produced following Sect 2 7 2 How to setup an HMGC run preparing the mode file TMODE and modi_inc The Fourier components considered in the simulation can be chosen considering the range of the q r profile typically starting from Mmin 1 up to Mmax ng r 1 For the specific case the modes m n 0 0 1 0 have been included for the equilibrium and the modes 1 2 up to 21 2 for the perturbation Mmax 2 x q r 1 2 x 10 5385 21 note that 4 X Mmax 1 80 is a v
29. 4 i 0 4 ue al 3 J gd x lt lt 6 Lj 45 Lu 45 LU Z Z z Lu Lu Lu 4 0 24 o 924 v IZ o Es Uv Uv a a a ag o o Q Q Q E 0 sl 0 Po E ZE 0 2 0 4 06 08 1 02 04 06 08 1 0 02 04 06 08 1 r a r a r a Figure 11 Frequency spectra during the linear growth phase twao 144 and dif ferent values for twindow first row twindow 48 wmin Wao 0 131 second row twindow 96 Wmin wao 0 0654 third row twindow 144 Wmin wao 0 0436 forth row twindow 192 Wmin wag 0 0327 ihann and ibuffer first column ihann 0 ibuffer 1 second column ihann 1 ibuffer 1 third column ihann 1 ibuffer 3 i e Wmin plot Wmin 3 49 f Gap Iplar Atog 1x 72 0 tow 432 00 Sum over m n l min l max 3 23 x S9AE US HMGC code ENEA Frascati 0 2 04 06 08 1 ria co f Gy Uo P At 1x 72 0 t7 432 00 Sum over m n min l max 3 23 X Q potisit B55E 05 0 24 0 14 HMGC code ENEA Frascat 0 2 0 4 06 08 1 r a Of Gap Iplar Atom 3x 72 0 t7 432 00 Sum over m n I min l max 3 23 389E 05 HMGC code ENEA Frascati 0 2 04 0 6 0 8 1 ria co f Gy Uo P At a 1x 144 0 top 432 00 w ogg plne Ato 1x 144 0 t 4352 00 co Gp o P At a 3x 144 0 tiop 432 00 Sum over m n I min l max 3 23 Sum over m n I min l max 3 23 Sum over m n I
30. HMGC User guide Last update June 13 2008 G Vlad S Briguglio G Fogaccia C Di Troia Associazione Euratom ENEA sulla Fusione C R Frascati C P 65 1 00044 Frascati Rome Italy Manuale utente di HMGC G Vlad S Briguglio G Fogaccia C Di Troia Riassunto Questo manuale descrive l utilizzo di HMGC Hybrid MHD Gyrokinetic Code il codice di simulazione 3 D ibrido magnetoidrodinamico girocinetico a particelle sviluppato a Frascati nei primi anni 90 HMGC stato sviluppato per studiare l interazione nonlineare di ioni energetici con turbolenza di tipo Alfv nico in plasmi che bruciano Il modello di plasma adottato nel codice HMGC consiste in una componente di plasma termico core e una popolazione di ioni energetici La prima descritta dalle equazioni della Magneto idro dinamica MHD ridotte O c nel limite di pressione nulla dove e l inverso del rapporto di aspetto del toro inclusi termini resistivi e viscosi La popolazione di ioni energetici descritta dall equazione di Vlasov girocinetica nonlineare espansa nel limite kipp 1 k essendo la componente perpendicolare al campo magnetico del vettore d onda e py il raggio di Larmor degli ioni energetici con gli effetti di orbita di deriva magnetica pienamente ritenuti e risolta con tecniche particle in cell PIC Lo scopo di questo manuale utente quello di rendere il lettore in grado di utilizzare il codice e di analizzarne i risultati con un insieme di st
31. Q 150 EPSILO 0 360781991d0 BUMPEQ 0 75D0 CG 0 2000D0 WG 0 220DO Fourier modes included in a HMGC simulation left the n 2 reference DIII D case right a n 2 4 case Black dots represent the modes actually included in the simulations the red crosses represent the modes considered in the simulation because of complex conjugate condition Radially integrated magnetic energy of the perturbed Fourier components vs time the flag fourth column in the TMODE_PL input file commands the Fourier components to be plotted Note that here the result of a simulation in which nout 100 is shown 6 4 2 alee di OI Example of integrated total energy Wiot m n vs time The labels on the right of the plot describing the m n poloidal and toroidal mode numbers of the curves are written every two modes because of space limitation although in the plot all the components are shown Radial profiles of the Fourier components from CTBO file as produced by the previous script only the first 8 frames are shown solid black line real part dotted red line imaginary part Example for density Gy red dashed curve is P g green dotted curve is die and ay radial profiles alu Bacal aia ine dai Frames number 161 twao 96 as produced by the program plot field f using the xplot field input data shown above Left 7 0 0 Centre r 0 y and superimposed the position of the first test particle prod
32. Sai 66 v SES x R 63 mH MH where I is the identity tensor I dij Fg t R M U is the energetic particle distribution function in gyrocenter coordinates and D 7 is the Jacobian of the transformation from 71 canonical to gyrocenter coordinates The Vlasov equation can be written as dFy 4 E o 64 with m d 0 d p di Ot dt gg and dZ dt given by the following equations of motion 14 20 21 22 dR Ub b b di Ub mes Vo LS x Vay E gue v dI 6x ving mg QH my dM 0 65 dt 65 dU La M b ae U allt V Va xVInB dt mg Og m MH QM a VajxVo I b5 vInB MHH MH Here ey and Qy eg B mgrc are respectively the energetic ion charge and Larmor fre quency The fluctuating potential a is related to the poloidal magnetic flux function through the relationship en Ro pr Note that the magnetic moment M is exactly conserved in this coordinate system and that ay 66 correspondingly neither Fy nor the equations of motion contain any dependence on the gy rophase 8 The particle simulation approach to the solution of Vlasov equation eq 64 consists in representing any phase space function G t Z by its discretized form G t Z PZG tZ 5 Z Z x DAG Z 6 Z Z 67 l where A is the volume element around the phase space marker Z and in assuming that each marker evolves in time according to the gyrocenter equations
33. Toroidal corrections enter the equations at the next order in the inverse aspect ratio expansion The derivation of these equations has been described in detail in ref 12 and is only briefly reported here Following the low tokamak ordering it is possible to write vi BI B B V N N go UA B V e po VY c VUE c ea VA VA vafa B at R Here a cylindrical coordinate system R Z y has been used and the subscript L denotes components perpendicular to Vy The magnetic field can be written as B F i F Vo RyV x Vy O amp B 50 where is the poloidal magnetic flux function Fp RoBo Bo is the vacuum toroidal magnetic field at R Ro and F O e Fp is given at the leading order by equilibrium 67 corrections Substituting eq 50 and Ohm s law eq 48 into Faraday s law eq 43 we obtain 1 Rot ve c m s B cV O ev4By en where is the scalar potential Taking the cross product by Vy eq 51 can be solved with respect to v cR Ro Bo Equation 52 states that at the lowest order the perpendicular velocity is given by the E x B Vox Vyt Oleva 52 vy drift Then taking the Vy component of eq 51 the following equation for the evolution of the poloidal magnetic flux function is obtained O cR c Og m Cung EE go 2 t AY O c vAB 53 with the Grad Shafranov operator A defined by eq 49 Upon applying the operator Vy
34. V x R to the momentum equation eq 40 the following equation for the evolution of the scalar potential is obtained a S 2c De Vier ve s c i ve Dt RoBo OZ Dt RoBo OZ Bo Bo 2 4 v4 By TBOVA Vo R VP e 4 To Y che x Ve O e 0 eni 54 where x R Porre oo pa DEA 10 9 Mes R 4 ROR OR 0O0Z Note that both in eq 53 and eq 54 v and F enter only at the fourth order in e In eq 54 the dependence on the density gradient has been retained explicitly With the particular choice of the mass density oR R const and using the definition of v given in eq 52 the continuity equation eq 39 is satisfied up to the third order The pressure equation becomes DP vA B o Dt Js 55 9 3 Hybrid MHD kinetic models 68 In order to include in this model the effects on an energetic ion population we can take advantage from the fact that the energetic particle density is typically much smaller than the bulk plasma density The following ordering can then be adopted 0 0 oe Heo where ny n and Ty T are the energetic particle bulk ion density and temperature respectively Thus the following ordering for the ratio of the energetic to bulk ion beta follows Bu Bi It can be shown 13 that making use of the above ordering the MHD momentum equation is modified by a term which represent the perpendicular component of the divergence of the energ
35. alid number for the ESSL FFT routine 7 3 How to setup an HMGC run plasma parameters file PARAM The bulk plasma density normalized profile should be assigned following one of the proposed expressions NPROFI flag see Table 3 For representing experimental data it is usually convenient to choose NPROFI 1 the three parameters a B pa can be obtained by fitting the experimental profile Normalized resistivity 7 and normalized viscosity v should be chosen such that numerical stability is assured typical values used in the simulations are 7 1 x 1075 v 1 x 1078 The initial amplitude perturbation AMP should be such to be well below the saturation amplitude in order to provide a sufficiently long linear phase of the simulation Actually if is initialized only a single Fourier component the one at the centre of the poloidal spectrum A time step DT must be chosen usually the MHD module loop on a shorter time step whereas the gyrokinetic module is called every NSUBCY time steps One should take care that the fastest particles considered in the simulations do not cross more than one cell in the 0 and directions or equivalent cell if the nogrid version is used to integrate correctly their motion 7 4 How to setup an HMGC run energetic particle profiles files den spli data and temp spli data The integration of the r q r profile will give a y r mesh see eq 9 in Sect 3 8 The Ecrit r Ecrito normalized profile
36. ational Atomic Energy Agency 8 Stix T H 1972 Plasma Phys 14 367 84 9 G Vlad F Zonca and S Briguglio Dynamics of Alfv n waves in tokamaks Rivista del Nuovo Cimento 22 7 1 97 1999 10 J P FREIDBERG Ideal Magnetohydrodynamics Plenum Press New York 1987 11 H R Strauss Phys Fluids 20 1977 1354 12 R Izzo D A MONTICELLO W PARK J MANICKAM H R STRAUSS R GRIMM and K MCGUIRE Phys Fluids 26 1983 2240 13 W Park S PARKER H BIGLARI M CHANCE L CHEN C Z CHENG T S HAHM W W LEE R KULSRUD D MONTICELLO L SUGIYAMA and R B WHITE Phys Fluids B 4 1992 2033 75 14 15 17 18 19 20 21 22 23 24 27 S Briguglio G Vlad F Zonca and C Kar Hybrid magnetohydrodynamic gyrokinetic simulation of toroidal Alfv n modes Physics of Plasmas 2 3711 3723 1995 G Vlad C Kar F Zonca and F Romanelli Linear and non linear dynamics of Alfv n eigenmodes in tokamaks Physics of Plasmas 2 418 441 1995 Also RT ERG FUS 94 16 Ass EUTATOM ENEA sulla fusione C R E Frascati Rome Italy 1994 G VLAD S BRIGUGLIO C KAR F ZONCA and F ROMANELLI in Theory of Fusion Plasmas edited by E SINDONI and J VACLAVIK pp 361 382 Bologna 1992 Editrice Compositori Societ Italiana di Fisica A BONDESON Nucl Fusion 26 1986 929 W W LEE Phys Fluids 26 1983 556 W W LEE J Comput Phys 72 1987 243 T S HAHM Phys Fluid
37. ax in u alpha if 0 0 use max available interval T iflag_contrai 1 sum from jerO_pl_min to jer0_pl_max 24 40 jerO_pl_min jer0_pl_max if not summed plot jer0_pl_min 0 iflag_relax 1 reset timkin_relax 120 timkin_relax_new 0 iflag_renorm 1 renormalize f_t_relax to n_H r t_relax and f t to n H r t only if iflag df 1 1 iflag filtro 1 plot of sedentary and displaced particles df only if iflag_df 1 5 51 er 0 375 0 625 tw y 0 00 a jer 0 375 0 625 twy 0 00 a jer 0 375 0 625 tw 0 00 j l l 9 uu fO x 5 812kb D 4 u f x 2 024 01 by x 6 417E 02 Ea t 3 1 u 3 S 1 2 n n n dA N N i i b o HMGC code ENEA Frascati 4 9 HMGC code ENEA Frascati b 9 HMGC code ENEA Frascati 1 o D o 0 5 0 5 1 Figure 13 Left fy r u u Centre fy r u u transformed on the E a mesh Right fu r E a All figures are obtained summing from rjero_pi_min 0 375 to rjero pl mar 0 625 Note that the input file xplot_deltaf_input shown in the text refers to the left plot the figure in the centre is obtained with iflag_versus 2 the figure on the right is obtained reading the file deltaf_ealfa data All figures refer to the frame 1 twao De u jer 0 375 0 625 t 47 126 00 a jer 0 375 0 625 tw y 126 00 a jer 0 375 0 625 twe 126 00
38. bution function The distribution function is stored on a r 14 u grid defined for plotting purposes More precisely the number of energetic particles Ngy ri hj x in the volume centered at r 4j uj is given by Nu ri nj ux np 0 8 fa ri uj ux dridujduy The file is in unformatted format The file deltaf ealfa data contains the same data but on a E o mesh instead of ju u The data stored on the ji u mesh can be transformed and plotted on a E a one note that this is not equivalent to store directly the distribution function in the E a space indeed the single particle energy dependence on 0 is lost The data can be also summed up between to radii to reduce particle noise see Fig 13 The plotting routines allow to subtract to the current fy r u u t the value of fp r p u t trelax where t trelax is a specific time thus the fx r u u is plotted see Fig 14 An example of the input file xplot deltaf input for the program plot deltaf f is listed hereafter deltaf data 0 ips 0 no PostScript file 1 PS file 2 EPS name follows 30 char pippo eps 1 321 1 ifirst step itot increm output time steps 1 iplO first plot 211 ipli last plot 1 iflag 1 plot deltaf 2 plot deltaf ealfa 0 iflag df 0 plot f 1 f t f t relax 1 iflag versus 1 mu u 2 E alpha only if iflag 1 0 00 0 00 am0 ami min max in mu E if 0 O use max available interval 0 00 0 00 u0 ul min m
39. cles in order to distribute their contribution to the pressure tensor components among the vertices of the cell Phase space coordinates and weights for the simulation particles are initially determined in such a way to yield a prescribed e g Maxwellian distribution function Particle pushing is performed by integrating eqs 65 by a second order Runge Kutta method more accurate than the standard O At Euler method O At is properly retained although more time consuming Particles that hit the wall r a are considered lost and are not re injected in the plasma 73 It has been shown 14 23 24 25 26 27 that as far as regimes are considered where the distribution function can be expected to slightly depart from the equilibrium one it is worth limiting the numerical investigation to the evolution of the perturbed part F p defined by the relationship Fy t R M U Fuo t R M U Fg t R M U 76 where F jo is the lowest order equilibrium distribution function In terms of F y eq 64 can be written in the form d F g r S 77 with Ts dme m Meanwhile eq 70 is replaced by the following one Npart D z 5Z Fu Z M wi t 6 Z Zi 0 78 1 with T t Aj6F pg t Zi t 79 and a wi AS t Zi t 80 Note that eq 76 is by no means equivalent to a linearization of the Vlasov equation since all non linear terms are correctly retained The decomposition of eq 76
40. ctric charge of the bulk ions plus impurity eventually and Te the electron temperature in MeV n and nj the electron and bulk ion densities Eg is the birth energy of the energetic particles and Ecrit o Ecrit 0 In the code the parallel velocity is normalized in this case to the energetic particle birth energy velocity u vino with vino VEo my 6 3 Bi Maxwellian distribution function idistr 3 The Bi Maxwellian energetic particle equilibrium distribution function fgQ i Maxw is defined as follows NHO mg 3 2 fEQ bi Maxw FEQbI Maxw 35 3 2 Ir T Ho an um exp MHU ks tcH 0 f amp Q bi Maxw non TL pry Q2 Tao 2T HoT y Timor v 56 Tru 4 Tiuori 37 Tn v Tinon 38 In this case the normalized quantities are u vino with ving VTLHo My and f uQeH0 T Ho The quantity oo is oo T1 go T go 57 7 How to setup an HMGC run In this section we will give indications on how to set up a specific run of HMGC we will refer to a DIII D discharge 122117 analyzed in Ref 7 W W Heidbrink provide us a set of radial profiles for the following quantities t 0 414 s of the DIII D discharge see Table 10 and Fig 17 ne r Electron density cm 3 ni r Thermal deuterium density cm 3 Nimp T Impurity density carbon cm 3 ny r Beam ion density cm 3 Note that this transp value is reduced by fast ion transport
41. dau damping on standard use of HMGC ICURV 0 curvature term off 1 curvature term on standard use of HMGC IOMST 0 w term off 1 w term on standard use of HMGC NPTEST number of test particles ITEST for each test particle enter 0 to initialize a true test particle 1 to follow a simulation particle 5 reals if ITEST 0 the initial coordinates t 0 of the test particle must or be given ee Brest Prost Htests Utest LTEST if ITEST 1 particle identification number ER_PERTO energetic particle pressure tensor term PREK set to zero for r lt ER PERTO standard use of HMGC ER_PERTO 0 1 iprek 1 PREK 1 0 l Fourier component inactive PREK 1 21 l Fourier component of the energetic particle pressure tensor term active defaults is all Fourier components active Table 8 Structure of the file KININP continued 2 29 3 8 Input files energetic particle density and temperature pro files Normalized to the on axis value energetic particle density profile and temperature if Maxwellian distribution function is loaded idistr 1 Eeri if slowing down distribution function is loaded idistr 2 or perpendicular and parallel temperatures if bi Maxwellian distribution function is loaded idistr 3 profiles must be provided on an equally spaced normalized poloidal flux function grid v These profiles are usually provided by standard transport codes e g available in the ITER database If experime
42. ear beam heated DIII D discharges between the energetic particle radial density profile expected from classical deposition and that deduced from the experimental measurements In spite of the slightly simplified physical model HMGC has been getting increasing attention from the international plasma physics community and it has been recently acquired by EPFL CRPP Lausanne University of California Irvine and IFTS Zhejiang University Aim of the present report is yielding a HMGC User Guide We proceed with a summary description of the various sections In Sect 2 it is shown how to produce a plasma equilibrium needed by HMGC In Sect 3 it is described the execution script which prepares the set of input files required for compilation and execution of the code In particular the script prepares both the sets of files required by the two modules that constitute the hybrid code the MHD module and the gyrokinetic one Sects 4 and 5 describe the output files of the MHD part and of the gyrokinetic one respectively they also describe the suite of graphics tools used to post process and visualize the results contained in these files Sect 6 describes the three types of energetic particles distribution functions that can be loaded to start a simulation the slowing down the maxwellian and the bi maxwellian distribution functions The various operations needed to setup a run of HMGC have been collected in Sect 7 where a specific HMGC run
43. ed to reduce e correction at the edge see Eq 8 namelist NRCHNL 6 diagnostic output channels at the radii RCHNL i 1 6 giving DIAPOS Real and Imaginary part of mn Tenni t Table 4 Structure of the file PARAM continued 19 The general time stepping in HMGC is as follow the MHD normalized time step is given by the parameter DT Each NSUBCY MHD time steps the gyrokinetic module is called This loop is performed NCYCLE times Thus a total umber of time steps of NTS NCYCLE NSUBCY is performed for each time sequence This time sequence is repeated NOUT times Thus the total number of time steps of a simulation is Ntime steps NTS NOUT NCYCLE NSUBCY NOUT and the total normalized time simulated is Ttotal DT NCYCLE NSUBCY NOUT Schematically using a fortran like schema these nested loops are as follows time 0 do i 1 NOUT do j 1 NCYCLE call GyroKinetic module do k 1 NSUBCY time time dt call MHD module enddo enddo enddo The time is normalized to the inverse of the on axis r 0 Alfv n frequency that is tcode t s wao s with wao vao Ro and vag the on axis Alfv n velocity If FREZOO true and FREZ10 true no evolution of equilibrium components and a sin gle perturbed n is included in the simulation a linear MHD simulation will be performed this is the usual operation condition of the MHD module The parameters epsil1 cgeps and cweps are used to define a radial function f r w
44. emperatures profiles e HMGC RESULTS directory containing output results 1 HMGC RESULTS caso equil DIII D 1 directory containing a DIII D equilibrium test case 2 HMGC RESULTS caso n2 DIII D 1 test 1 directory containing a DIILD test case 65 9 Generalities on HMGC Beginning of excerpt form the Appendices of Ref 9 For the references to equations not resolved please refer to the original paper 9 1 MHD equations We wish to start the discussion of the general properties of the Alfv n wave spectrum using simple model equations the so called ideal magnetohydrodynamic MHD equations which in Gaussian units read 10 Jo 90 v 39 at t V v 0 39 dv 1 4 27 VP JxB 40 d P Xe M aes 41 zr S an bd ov lt BR L0 42 C 10B Ec ccr 4 V x dec 43 Veg zy 44 C V B 0 45 In the above equations v is the fluid velocity J is the plasma current B is the magnetic field o is the mass density P is the scalar pressure of the plasma T is the ratio of the specific heats c is the speed of light and d MES 4 di ob V 46 is the convective derivative The ideal MHD equations describe the plasma as a single fluid In particular eq 39 describes the time evolution of mass conservation of the total number of particles Equa tion 40 describes the time evolution of momentum showing that the fluid is subject to inertial pressure gradient and magnetic forces Equation 41
45. etic particle stress tensor Iy in ref 13 an alternative equivalent form in which the electric current associated to the energetic ions appears instead of the energetic particle stress tensor is also derived Thus the O c equation for the evolution of the scalar potential is modified as follows D 2c Od 9 i D c Od B 5 zm Vie vo 5 zz V aS pega Oy PANI x Ve e viBe 56 Arce cRo gt d d ale In order to close the set of reduced MHD egs 53 and 56 the hot particle stress tensor components can be evaluated by directly calculating the appropriate velocity moment of the distribution function for the particle population moving in the perturbed fields v and see appendix sect 9 4 9 4 Hybrid MHD kinetic code In this section we describe the code that solves the O c reduced MHD model in the limit of zero bulk plasma pressure In such limit only eqs 53 and 56 need to be solved As a boundary condition we take a rigid conducting wall at the plasma edge The numerical tool 14 15 16 used to solve the O e model is based on a field solver originating from an existing O e code 17 Such field solver uses toroidal coordinates r v finite differences in the radial direction r and Fourier expansion in the poloidal 9 and toroidal vy direc tions The coupled equations for the Fourier components of the magnetic and velocity stream functions w and are advanced in time using a semi
46. ferent values of MODECH will provide plots of other quantities and will require different input data 33 eJaok031 CIIO DII D caso sldown n2 Dill D 1 06 33 Hol bia nem NIN 3 CONMORW PCOWODMIMOARWNY 0X INS IND NI NI INS NO NI NA NO NA NA RS SO SONO SONORO NO n nus ggg ISS or NU Figure 4 Radially integrated magnetic energy of the perturbed Fourier components vs time the flag fourth column in the TMODE PL input file commands the Fourier components to be plotted Note that here the result of a simulation in which nout 100 is shown 34 After running the program epe3ak31 u an output file is written which can be used to produce a time sequence of plots movie by the fortran plot_energy f An example of the input file xplot_energy_input for the program plot_energy f is listed hereafter ENERGY TMODE_PL 0 ips 0 no PostScript file 1 PS file 2 EPS name follows 30 char pippo eps 15 10 ifirst step itot increm output time steps 3 i magnetic energy 2 kinetic energy 3 magnetic kinetic energy 0 time_max if 0 take time_max from stop_run file 1 e 32 0 y_en_min y_en_max 161 iplO first plot 161 ipli last plot 5 The above input data example contains the name of the input data files ENERGY TMODE PL ips is a flag for producing a PostScript output followed by the name of such a file pippo eps ifirst step itot are the first and last o
47. h and possibly destabilize Alfv n modes Energetic ion transport and confinement properties of crucial importance for achieving efficient plasma heating and therefore ignition conditions can in turn be affected by nonlinear interactions with the Alfv nic turbulence Thus large efforts have been devoted to assess the stability of shear Alfv n modes in tokamaks and to investigate their effect on the energetic ion transport The need of fully retaining nonlinear dynamics and properly taking into account kinetic effects such as resonant interactions between energetic ions and Alfv n modes and the nonperturbative character of such interactions makes the numerical particle simulation approach the natural tool for this investigation A hybrid MHD particle 3 D simulation code HMGC Hybrid MHD Gyrokinetic Code has been developed in Frascati in the early 90 s The plasma model adopted in the HMGC code consists of a thermal core plasma and an energetic ion population The former is described by reduced O e Magneto Hydro Dynamics MHD equations in the limit of zero pressure e being the inverse aspect ratio of the torus including resistivity and viscosity terms The reduced MHD equations expanded to O e allow us to investigate equilibria with shifted circular mag netic surfaces The energetic ion population is described by the nonlinear gyrokinetic Vlasov equation 1 2 expanded in the limit ki px 1 with k being the perpendicu
48. hich modulates the toroidal corrections ii de 22 ftam ie T l fer 1 fo r lt ceg 8 Ro fe for r gt Cgc 20 3 3 Input file stop run At the beginning this file contains only one line which can be used to overwrite the parameter NOUT given in the file PARAM 20 nout_new At run time the same file is written and read by HMGC and this allows to stop or extend the run by editing it and changing the parameter nout 25 nout 21 ncount 192 0000000 time Here ncount is the actual number of time sequences performed by the code and time is the corresponding normalized simulation time 2l 3 4 Input file TMODE This input file contains the list of the Fourier modes actually included in the simulation The first column is the mode index running from 1 to LM In the second and third columns are listed the corresponding poloidal m m and toroidal n n I mode numbers It is assumed that the modes are ordered by increasing n and for each n by increasing m ars COONDORWNE IP BR FOWOONDOPWNFEFO 12 NRPRPRPRPRPRPREH COOANDOORWND NORPRPRREPREFR Od 00 10 01450 NIN Ne NNNNNNNNNNNNNNNNNNNNNCO N w N arg 22 3 5 Input file TMODE PL This input file is used by the plot post processing program epe3ak31 u It is exactly the same than the input file TMODE with the add of a fourth column in which 0 means that this component will
49. ical energy for slowing down distribution function IDISTR 2 ALF parameter for non uniform energetic particles radial loading uni form fraction 0 lt ALF lt 1 NPIC parameter for non uniform energetic particles radial loading num ERO I DEL I ber of gaussians overimposed to the uniform distribution fraction ALF If NPIC ne 0 the corresponding values of the radial positions and widths ERO I DEL I of the gaussians must be given Table 6 Structure of the file KININP 27 Quantities Comments ENHSNI nyo N o on axis ratio between the energetic particle and bulk ion densities EMHSMI mg m ratio between the energetic particle and bulk ion mass timkin_anu parameter for ramping the energetic particle density energetic par ticle density ramps for t gt timkin_anu timkin_anu should be greater than TIMKIN RELAX but smaller than the time at which non linear phase occurs Set timkin_anu greater than the total simula tion time to avoid ramping TIMKIN_RELAX time at which the energetic particle distribution function is assumed to be relaxed because of initialization of the distribution function in terms of non conserved quantities ANU_MAX parameter for ramping the energetic particles multiplying factor of the normalized energetic particle density ANU_DOT_O parameter for ramping the energetic particles time derivative of the normalized energetic particle density i_write_deltaf
50. implicit algorithm where all the linear terms that couple with the cylindrical part of the equilibrium i e the component having 69 poloidal and toroidal mode numbers m n 0 0 are treated implicitly The non linear terms the terms which arise from the toroidal corrections to the cylindrical approximation and the contributions of the energetic particles the term containing V IIg in eq 56 are treated explicitly Moreover only the Fourier components in a half plane of the m n space are evolved the ones that fall in the other half plane being recovered from the reality condition of the solution Wy alts t docu t EE t Peni t 57 The equilibrium configuration used for numerical simulations can be exactly calculated to the desired order in e starting from the expression for the equilibrium toroidal current J 58 c Ro Op and expanding v4 in powers of e P r 9 po r Wi r 9 O v In the toroidal coordinate system r V p the Grad Shafranov operator can be expressed as A18 lo 1898 1 g2 sind o rOr V ar 1299 RV 8r r Ov with R Ro r cos 9 To the leading order that is in the cylindrical approximation eq 58 is given by ld dig N 4m E yielding dy ES Bo r dr Roqlr which can be integrated assigning Y a 0 and g r the safety factor in the cylindrical O e approximation g r rBog RoBoy Equation 59 gives the symmetric m 0 n 0 Fourier co
51. ized resistivity i e the inverse of the Lundquist number S the ratio between resistive and Alfv n times S r 74o with 7 uoa n and Tao w49 VISC similar to ETA parameter but representing viscosity DT elementary time step NPRI some outputs are performed every NPRI time steps NPRI must satisfy NTS NPRI lt MAXPRI RWALL resistive wall normalized radius parameter required but not used by fortran give any real number TAUWAL resistive wall characteristic time parameter required but not used by fortran give any real number VEDGE plasma bulk velocity at the edge parameter required but not used by fortran give any real number CURAMP current ramp only significant if equilibrium is evolved FREZOO logical variable if true does not evolve freeze 1 1 mode 0 0 FREZ10 logical variable if true does not evolve freeze 1 2 mode 1 0 EQUIL logical variable if true the code is used to compute an equilibrium useful for nonlinear MHD runs DROP logical variable if EQUIL true and DROP true kinetic energy is removed by dropping m n by a fixed factor useful for preparing nonlinear MHD runs Table 3 Structure of the file PARAM 18 Quantities Comments NPROFI switch for assigning the bulk density radial profile A NPROFI 0 1 NPROFI 1 f pla b pa NPROFI 2 p g r g 0 aligned toroidal gap ALFA
52. li lit o 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 0 0 1 0 2 0 5 0 4 0 5 0 6 0 7 0 8 6 8 1 Figure 1 ITER SC4 q profile example of use of current bumps left no bumps right with two bumps one bump is positive in amplitude BUMPEQ 1 30D0 located at r CG 0 4000D0 r 0 632 having width WG 0 350DO and the second is negative bumpeqi 2 20D0 at r cgi 0 90D0 r 0 949 having width wg1 0 30D0 In figure 2 the q profile used to simulate the DIII D discharge 122117 at t 0 414 s is shown Hereafter it follows the EQUIPA namelist used amp EQUIPA Qo Qi 3 9891D0 15 000DO 11 RL 2 5D0 NREQ 150 NMESHA 0 NPOIDA 2 APLACE 1 0 426D0 0 9DO 0 00DO 0 0 0 0 0 0 0 AWIDTH 1 0 100D0 0 10D0 0 00DO 0 0 0 0 0 0 0 SOLPDA 0 60D0 EPSILO 0 360781991d0 RHOFLG FALSE BETAO 0 00000D 0 C1 1 7438D0 C2 2 3515D0 C3 12 01D0 C4 15 988DO C5 7 3964DO BUMPEQ 0 75DO CG 0 2000DO WG 0 220DO bumpegi 0 00DO cgi 0 90D0 ugi 0 30DO bumpeg2 0 00DO cg2 0 95DO wg2 0 20d0 ireadcur 0 amp END ITER DIII D 1 10 a L B amp L L 4 L o amp IE e uL 0 birt lle bret 0 0 1 02 0 5 0 4 05 06 07 0 8 0 9 1 Figure 2 DI D discharge 122117 at t 0 44 s The parame ter used are QO 3 9891D0 Q1 15 000D0 RL 2 5D0 NREQ 150 EPSILO 0 360781991d0 BUMPEQ 0 75D0 CG 0 2000D0 WG 0 220DO
53. mic MHD Alfv n waves Tokamaks Particle in cell PIC techniques Gyrokinetic simulations Contents 1 Introduction 2 How to produce an MHD equilibrium file 3 Execution script 3 1 3 2 3 3 3 4 3 5 3 6 3 7 3 8 Include file modi is 254 2 3 1 o Rete od e os ed Lee Input file PARAM eee A i vod Gee ht ib a Inputs tle stop rum s su seta busa barre ibra aerea Ge a Input file TMODE lire de pala ors pct ee ee ee ee e Tnp t file TMUDESPL ac alata 2 4 adatta of aida d api Include Mle RIA die enea phon a te d IR dene d IR det d Rut Input file KININP o ut Rue rhet ace cjr roe Seco Ree NEC Input files energetic particle density and temperature profiles 4 Output files produced by the MHD module 4 1 4 2 4 3 Outp tfMe CTO csi dela d ee gba Bag usse qos elc x A Output file ENERGY cc toon Gob Mien P bonded hd eines d edP ae pes Output nie TBU eus nieder en de dece Se o e dec Pret de oen de e es 5 Output files produced by the gyrokinetic module 5 1 5 2 5 3 5 4 5 9 Output file OUTDNT riore tee tek die Output file TESTWRIK s parearea ae RE ele i ea Outputs file PHIWRITE and APWRITE 04 Output file deltaf data deltaf ealfa data Output file power data ue an oe uo eek ORS OO RR SES 6 Energetic particle distribution functions 6 1 6 2 6 3 Maxwellian distribution function idistr 1 Slowing down distribution function idistr 2 Bi Maxwellia
54. min l max 3 23 i Xi S85E 05 E ret XS CZIE GS ri X TAM00E 05 A x 9 UV a A mn 2 2 0 24 T 0 24 T 0 24 T E d i z Lu Lu Lu 0 14 Lot Lot 4 Es La a a a a o o Q Q Oo Oo i 0 i o i 0 2 04 06 08 1 0 02 04 0 6 0 8 1 0 02 04 0 6 08 1 rf r rf co f cog o P Ato 1x 288 0 ti 432 00 co plp Ato 1x 288 0 tw 432 00 co f cog Go P Aton 3x 288 0 t7 432 00 Sum over m n I min l max 3 23 Sum over m n I min l max 3 23 Sum over m n l min l max 3 23 et ra X 8B 9E US eX 7 BBEROS iii X 9 34 BE UD HMGC code ENEA Frascati 1 0 2 0 4 0 6 0 8 r a HMGC code ENEA Frascat 0 2 04 06 08 1 r a HMGC code ENEA Frascati 0 2 04 0 6 0 8 1 r a Figure 12 Frequency spectra during the saturated phase tw4 432 and differ ent values for twindow first row twindow 72 Wmin Wao 0 0873 second row twindow 144 Wmin wao 0 0436 third row twindow 288 Wmin wao 0 0218 ihann and ibuffer first column ihann 0 ibuffer 1 second column ihann 1 ibuffer 1 third column ihann 1 ibuffer 3 i e Wmin plot Wmin 3 50 5 4 Output file deltaf data deltaf ealfa data The file deltaf data contains the time series of the energetic particle distribution function fu r p u thus integrated over 0 and is read by the program plot deltaf f see Sect 6 for the details on the distri
55. module kin uncomplete F GK module extr push complete F pushing routine of GK module extr pressure complete F pressure routine of GK module commons for HMGC commr31 input commr31 uncomplete F double calcolo nne f program to initialize the particles in the E a space interp spline f spline program to interpolate experimental data on the HMGC mesh nlr_interp_spline common for the spline program interp spline f psi from rho q exp f program to generate r from r q r experimental data upda f simple pre processing program to produce goofy f programs from goofy u upda ksh shell script to run upda pwr5 version temporary directory architecture dependent for producing the executable makefile makefile 3 HMGC graphics directory containing graphics post processor programs a b c d NM SAA VALI epe3ak31 f profilk f plot density f plot energy f 64 plot field f plot deltaf f g plot power f plotta_max1 f plot equil f simple plotting program to compare a computed equilibrium with experimental data j GVGRAPHOLIB directory containing some graphical routines library libgvgraph0 a k GVGRAPHILIB directory containing some graphical routines library libgvgraphi a e HMGC INPUTS 1 HMGC INPUTS equilibrium directory containing equilibrium files 2 HMGC INPUTS profiles directory containing energetic particle experimental density and t
56. mponent of the poloidal magnetic flux function v To the next order in e eq 58 yields r Or Or coy 10 Ovi 10242 1 dig _ r2 092 Ro dr 4T d R vet J 60 C 1 dye4 Ro De yea po Equation 60 admits solutions of the form eq EMT di 1 1 0 v4 r cost A r cos qu 61 70 where we have introduced the so called Shafranov shift A Substituting eq 61 into eq 60 and using the leading order solution of eq 59 the following equation for the Shafranov shift A is obtained A ld dig dA 1 de 0 2 r dr 87 dr x 2 i 62 Equation 62 can be integrated assigning the radial derivative of the Shafranov shift at the center A 0 0 regularity condition and A a 0 corresponding to v 0 on the rigid conducting wall to obtain A r The substitution of A r into eq 61 allows us to obtain the first order 1 0 Fourier component of the magnetic flux function v thus completing the equilibrium solution at the desired order Note that once fixed r and 0 the quantity A r corresponds to the shift with respect to the center of the poloidal cross section of the geometric center of the magnetic surface labelled by the value v r 9 Such a shift causes shear Alfv n waves even when propagating along the magnetic field line to cross the radial grid thus imposing restrictions on the time step of integration 12 Further restrictions are imposed by the strength of the ex
57. n distribution function idistr 3 7 How to setup an HMGC run Tal How to setup an HMGC run preparing the equilibrium file EQNEW 9 15 16 17 21 22 23 24 26 30 31 31 31 37 40 40 43 45 51 53 54 54 56 56 58 61 7 2 How to setup an HMGC run preparing the mode file TMODE and modi inc 61 7 3 How to setup an HMGC run plasma parameters file PARAM 61 7 4 Howto setup an HMGC run energetic particle profiles files den spli data and temp spli data bat Salito raid I 61 7 5 How to setup an HMGC run energetic particle dimensioning file grid_inc 62 7 6 How to setup an HMGC run energetic particle parameters file KININP 62 HMGC directories structure 64 Generalities on HMGC 66 9 1 lt MED equationz deo h Hess died cede ido e coe 66 92 Order e reduced MHD 5300 08 dela aldo 67 9 8 Hybrid MHD kinetic models 1 8 Duae eS ae tea La eta 68 9 4 Hybrid MHD kinetic Coder ue uou a ees Soe ee es wie du 69 List of Figures 1 ITER SC4 q profile example of use of current bumps left no bumps right with two bumps one bump is positive in amplitude BUMPEQ 1 30D0 located at r CG 0 4000D0 r 0 632 having width WG 0 350D0 and the second is negative bumpeq1 2 20D0 at r cg1 0 90DO r 0 949 having width wgi 0 3000 DIII D discharge 122117 at t 0 414 s The parameter used are QO 3 9891D0 Q1 15 000D0 RL 2 5DO NRE
58. n specie D Ay 2 Zp 1 toroidal field Bro 2 T Ro 1 688 m a 0 609 m Carbon impurity Ey 0 077 MeV volume averaged value of ppar Pperp 1 44 for the beam ions 62 RHOSA fana 0 032863457d0 1 02 x 10 L AY Zp Ey ev Brest n Vin 095010 Bo solae U AO 2 18 x 105 Bro rea X Aini 1020xm 3 Asbi VTHSVA Ui A0 0 271063836d0 4 491 Ej s y mg Boss X Aini 1020 xm 3 ale A imp mp 1022xm 3 usdelta_0_input 1 A 2 3256d0 The corresponding value of A has been chosen in order to match the experi mental value of the ratio ppar Pperp 1 44 cosalfa_0 input COS Qo 0 68128d0 Rtan Ro 1 15m 1 688m e0sec0 Eo Ecrito 4 153850158d0 Ecrito from Eq 34 evaluated at r 0 ENHSNI nuo Nio 0 264848976d0 ngo n Nimp EMHSMI mg mi 1 d0 IDISTR 2 slowing down distribution function loaded Table 12 Preparing the KININP file 63 8 HMGC directories structure e 7HMGC main root directory 1 HMGC script directory containing compilation and execution scripts xe3_HMGC execution script of e3_complete F xequil_HMGC execution script of eqe3aaab u xplot HMGC execution script for post processing plot programs e g epe3ak31 profilk plot density plot field plot deltaf plot power xpsi from rho q exp execution script of psi from rho q exp f 2 HMGC sources directory containing HMGC sources n eqe3aaab u equilibrium module e3 complete F MHD
59. normalizations in the gyrokinetic module will be different The expressions shown in egs 1 2 3 4 are appropriate for describing a mono tonic q profile but they are inadequate to describe more general safety factor profiles as e g reversed shear profiles Thus a number of bumps on the current density profile can be superimposed on eq 3 Actually up to 3 bumps can be superimposed 5 r2 cgi i l1 f Via vn x 2A X do s 5 where bump can be positive or negative The current density profile resulting from Eq 5 is then rescaled and such to provide a q profile which has the minimum equal to the parameter qo of Eq 1 The meaning of the different parameters of the input file EQUIPA assigned as a namelist with the same name of the input file are listed in Table 1 QO minimum q value Qi maximum value of q at r a 1 qa if bumpegi 0 RL NREQ Number of points in the radial mesh NMESHA parameter of non equally spaced mesh usually not used NPOIDA parameter of non equally spaced mesh usually not used SOLPDA parameter of non equally spaced mesh usually not used APLACE i parameters of non equally spaced mesh usually not used AWIDTH i parameters of non equally spaced mesh usually not used EPSILO inverse aspect ratio e a Ro RHOFLG logical value if true compute n r mo jo j r BETAO parameter for equilibrium pressure pr
60. not be considered by the plot program whereas 1 means that it will be considered 1 0 0 0 2 1 0 0 3 1 2 1 4 2 2 1 5 3 2 1 6 4 2 1 T 5 2 1 8 6 2 1 9 7 2 1 10 8 2 1 11 9 2 1 1210 2 1 1311 2 1 14 12 2 1 15 13 2 1 16 14 2 1 17 15 2 1 18 16 2 1 1917 2 1 2018 2 1 2119 2 1 2220 2 1 2321 2 1 23 3 6 Include file grid inc This lines of fortran are included in the MHD module e3_complete F and in the gyrokinetic module trough the general common commr31 complete F of HMGC It defines general particle simulation parameters see also Table 5 Note that some of the following parameters refer to toroidal and poloidal meshes also if the code is run in the nogrid mode see Ref 6 Particles here are simulation particles The file grid inc is constructed from the two files nlr_inc and altri _grid_inc_1 written in the execution script cat gt HOMEroot_sources pwr _version nlr_inc EOF PARAMETER n1r 64 EOF cat gt HOMEroot_sources pwr _version altri_grid_inc_1 lt lt EOF PARAMETER NTH 168 amp nintphi 2 nmax 1 amp nph_su_nintphi 4 amp NPH nph_su_nintphi nintphi amp ippc 2 amp nne 672 amp npart nlr nth nph ippc 3 amp NMODOM 27 amp NRZ 5 EOF Please note that NLR should be NLR lt NR see Sect 3 1 where NR is defined NTH should be chosen such that NTH gt 2Mmax see Sect 3 in the following example NTH 8Mmax
61. ntal profiles are provided e g p q p nn p Telp Tu p as functions of p Pjimiter the usual radial like co ordinate of transport codes e g TRANSP with the toroidal magnetic flux function a simple program is provided psi from rho q exp f which returns the poloidal flux function yw in terms of p integrating the following expression _ d m alp to obtain v v p The normalized coordinate proportional to the poloidal flux function should be such that it is zero in the centre and unity at the edge Those profiles will be interpolated using splines on the desired equally spaced normal ized v mesh by the fortran program interp spline f Experimental files with Wnorm p NH pnorm P and similar for the other quantities must be provided their names are e g den exp DIII D 1 temp exp DIII D 1 temp exp par xxx where temp exp DIII D 1 contains the energetic particle isotropic temperature in the case idistr 1 the Eeri normalized profile in the case idistr 2 or the energetic particle perpendicular temp exp DIII D 1 and parallel temp exp par xxx temperature profiles in the case idistr 3 Then a corresponding output on the equally spaced norm mesh will be produced 30 4 Output files produced by the MHD module 4 1 Output file CTTO The CTTO file has exactly the same form of the INCOND file but includes all the Fourier com ponents 4 2 Output file ENERGY The ENERGY file contains the fo
62. o 1m do not enter min max 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT 38 0000 r o E S p 0000 r o 0 6 0 8 1 960 0000 r o t 960 0000 r o 0 2 04 0 Figure 6 Radial profiles of the Fourier components from CTBO file as produced by the previous script only the first 8 frames are shown solid black line real part dotted red line imaginary part 39 5 Output files produced by the gyrokinetic module 5 1 Output file OUTDNT This file contains the time series of the radial profiles of energetic particle density DNTOT nor malized energetic particle density perpendicular PPERP and parallel PPARA energetic particle pressure The gyrokinetic module produces an output on the file OUTDNT every NWRITE time steps in analogy with the MHD module with NWRITE taking the place of NPRI thus the output time step is given by Atoutput Gg DT NWRITE in Alfv n time units WRITE 43 ISTEPO TIMKIN DENOUT IF ISTEPO EQ 0 THEN WRITE 43 ASPECT WRITE 43 USPS00 write 43 0 WRITE 43 NERRE LMEQ DO IR 0 NERRE WRITE 43 RLOWCIR ENDDO DO IR 0 NERRE DO LMODE 1 LMEQ WRITE 43 PSIEQ LMODE IR ENDDO ENDDO DO LMODE 1 LMEQ WRITE 43 MMODE LMODE ENDDO ENDIF DO 1 JJER 0 NLR WRITE 43 DNTOT JJER ANTOT PPERP JJER PPARA JJER 1 CONTINUE ISTEPO is the current time step index TIMKIN is the normalized current time DENOUT 0 0 USPS00 wW0 0 7 ma
63. of zeros can be added to the time sequence using the parameter ibuffer Note also that the plotting routines interpolate the resulting FFT function in order to draw the contour plot curves care should be taken in choosing the twindow parameter in order to resolve a specific mode see Fig 11 and Fig 12 to compare the effect of different parameters on the produced spectra ey 0 000 top 96 00 oxy 5 306 t 57 96 00 twy 96 00 x amp D44E D09 testparticle 1 zp4gp pg testparticle 1 3 9 1 3 os 9 os 9 n go o5 LL o5 LL S LL sa ga ta 5 Ed lt V lt 2 Lu amp 2 Lu t uJ o E o a D lt a2 5 a V o g di c ES 8 Qu i 0 as O as 0 n a ET 00 20 5 Q am 43 0 o gt I TE I taa Figure 8 Frames number 161 tw49 96 as produced by the program plot field f using the xplot field input data shown above Left r 0 0 Centre r 0 and superimposed the position of the first test particle produced by assign ing iflag test 1 Right trajectories of the first test particle rist test tests Utest produced by assigning iflag_test 1 and iflag_trajectory 1 47 ton 96 00 tw p 96 00 test particle 1 poloidal plane test particle 1 equatorial plane HMGC code ENEA Frascati HMGC code ENEA Frasca Figure 9 Test particle trajectory projected on the poloidal cross section y 0 the cross indicates r 0 and on the equatorial plane
64. ofile usually not used C1 C2 C3 C4 C5 parameters for equilibrium pressure profile usually not used BUMPEQ CG WG bumpeg 1 cg1 wgi BUMPEQ1 CG1 WG1 bumpeq C92 wg BUMPEQ2 CG2 WG2 bumpeg 3 C93 W93 ireadcur parameter to read current density profile as alternate input Table 1 Parameters in the file EQUIPA 10 To help in fitting an experimental q profile an utility to compare the experimental profile with the one obtained with the program eqe3aaab is provided program plot_equil see fig 1 In figure 1 the effect of including or not including the bumps in the current profile is shown amp EQUIPA QO 2 4110D0 Q1 5 1280D0 RL 4 0D0 EPSILO 0 293217665d0 BUMPEQ 1 30D0 CG 0 4000D0 WG 0 350D0 bumpeqi 2 20D0 cgi 0 90D0 wgi 0 30D0 amp END Please note that peculiar shaping of the current density profile should be avoided as much as possible in order to prevent the not desirable growth of MHD unstable modes Note that a positive bump in the current profile is used to produce an off axis minimum in the q profile whereas a positive bump at the edge is used to pull up the q profile at the edge ITER SC4 no bumps ITER SC4 bumps E E sl L amp E L s si E B ak L 4 E k L E E L L 3L C 3 LE E r q E 2L IC 2 Sib L E E C 1 F ru l fosse left L L j UT t E Pan NEU DI Lae L o FI A E A ETTI TOI Cri ilivialiv ii Pei ii dai via
65. or the program plot power f is listed hereafter ENERGY power data 0 ips 0 no PostScript file 1 PS file 2 EPS name follows 30 char pippo eps 1 321 10 ifirst_step itot increm output time steps 161 iplO first plot 161 ipli last plot 2 iflag_versus 1 mu u 2 E alpha 0 00 0 00 am0 ami min max in mu E if 0 0 use max available interval 0 00 0 00 u0 ui min max in u alpha if 0 0 use max available interval 1 iflag_contrai 1 sum from jerO pl min to jer0_pl_max 24 40 jerO_pl_min jerO pl max if not summed plot jerO pl min 5 u jer 0 375 0 625 twe 96 00 jer 0 375 0 625 twy 96 00 uu power x 2 159E 14 m u power x 2 954E 14 34 l 4 OB 2 a N pole 5 5 T 4 o HMGC code ENEA Frascati HMGC code ENEA Frascati o l Figure 15 Left P r u u Also the curves representing the maximum energy loaded in the initial distribution function dotted lines and the trapped untrapped particle boundaries for the lower solid line and upper dashed line radii considered are shown Right Pg r E o All figures are obtained summing from fjero pimin 0 375 to Tjer0_pl_max 0 625 53 6 Energetic particle distribution functions HMGC considers at present three different initial energetic particle equilibrium distribution func tions fgg namely anisotropic Maxwellian idistr
66. plicitly solved terms as e g in the case of high inverse aspect ratio equilibria and or highly non linear cases The term IIy in eq 56 is the pressure tensor of the energetic hot ions it can be expressed in terms of the corresponding distribution function fy Ig my f d vvv fy with my being the energetic ion mass to be determined by solving the Vlasov equation the collisionless limit of the Boltzmann equation Since the time scale of the dynamics we want to analyze is long compared to a cyclotron period it is convenient 18 19 to solve the Vlasov equation in the gyrocenter coordinate system Z R M U 9 where R is the gyrocenter position M is the magnetic moment U is the parallel velocity i e the velocity along the magnetic field line and is the gyrophase This corresponds to averaging the single particle equations of motion over many cyclotron orbits and allows one to retain the relevant finite Larmor radius effects without resolving the details of the gyromotion Such a choice is particularly suited for numerical time integration of the particle motion as the numerical stability constraint on the time step size turns out to be much less severe than that we would obtain without adopting the averaging procedure The hot particle pressure tensor assumes the following form in terms of the gyrocenter coordinate distribution function Ze 1 _ Li cli Ily t x 7 d6ZD _zFu t R M U x MH QyM_ a 2 QyM
67. rapped particle boundaries for the lower solid line and upper dashed line radii considered are shown Right Py r E a All figures are ob tained summing from Fjero_pimin 0 375 tO Fiero pt moz 0 625 Pri DA hws ou ddl adr do d T S oe erariale Experimental profiles for the DIIT D shot 122117 Beat geometr y se go So Ys eet ea eni er e ae A a List of Tables Oo ON So om 4 W SO I e e Ne c Parameters in the file EQUIPA 2 000084 10 Structure of the file EQNEW LL 13 Structure of the file PARAM 22 5 no Rem Phe ERE a 18 Structure of the file PARAM continued stopper oe lele 19 Structure of the file grid_inc esas ek eth ee tre pioli ae DER ee 25 Structure of the file KININP s Sor xk ER oe hes BS Es Re D ES 27 Structure of the file KININP continued 4 iu ridi Xd 28 Structure of the file KININP continued 2 29 Quantities in the file ENERGY x vac baci ae ea 32 Experimental radial profiles provided by DI D team 58 Preparing the grid_inc file scafo grato ne dela iena di 62 Preparing the KININP file 4 225 xe et og Rose EE ws 63 HMGC User guide 1 Introduction In a burning plasma fast ions MeV energies produced by additional heating and or fusion reactions are expected to transfer their energy via Coulomb collisions to the thermal plasma particles 10keV energies Due to their high velocity of the order of Alfv n velocity they can resonate wit
68. referred to a DIII D beam heated discharge is used as an example Sect 8 shows the list of the directories tree structure of HMGC Finally in Sect 9 several excerpts of Ref 9 are reported to illustrate the analytical details of the model that constitutes the basis of HMGC 2 How to produce an MHD equilibrium file The equilibrium file required by HMGC is produced by running the fortran file eqe3aaab This program solves the Grad Shafranov equation expanded to the O e in the inverse aspect ratio e a Ro with a and Ro the minor and major radius of the torus respec tively for the poloidal flux function w see Sect 9 assuming an analytic parametriza tion of the safety factor profile q q r with r the normalized minor radius r r a a being the minor radius of the circular cross section torus given by 1 2 n l 1 with ro defined in terms of A and the q value at the centre q r 0 qo and at the q r 4o edge q r 1 qu 1 2A ees 8 2 The normalized to Br Ro current density profile and the shear profile can be derived from the previous expressions je i 8 Eo i TO From the above expressions the normalized to B7a Ro Fourier components 1 o for the equilibrium poloidal flux function are obtained namely wo V here m is the poloidal mode number and the toroidal mode number n 0 has been assumed 9 for the axisymmetric equilibrium Please note that the
69. rumenti grafici anch essi descritti con un certo dettaglio Parole chiave Prodotti di fusione Particelle alfa particelle veloci Magnetoidrodinamica MHD Onde di Alfv n Tokamaks Tecniche Particle in cell PIC Simulazioni girocinetiche HMGC user manual Abstract This user guide describes the use of HMGC Hybrid MHD Gyrokinetic Code the hybrid MHD particle 3 D simulation code developed in Frascati in the early 90s HMGC has been written in order to study nonlinear interactions of energetic ions with the Alfv nic turbulence in burning plasma conditions The plasma model adopted in the HMGC code consists of a thermal core plasma and an energetic ion population The former is described by reduced O c Magneto Hydro Dynamics MHD equations in the limit of zero pressure e being the inverse aspect ratio of the torus including resistivity and viscosity terms The energetic ion population is described by the nonlinear gyrokinetic Vlasov equation expanded in the limit ky pu lt 1 with k being the perpendicular component of the wave vector to the magnetic field and pg the energetic ion Larmor radius though fully retaining magnetic drift orbit widths and solved by particle in cell PIC techniques The aim of this user guide is to make the reader able to run the code and analyze its results using a suite of graphics tools also described in some detail Keywords Fusion products Alpha particles Fast particles Magnetohydrodyna
70. s 31 1988 2670 T S HAHM W W LEE and A BRIZARD Phys Fluids 31 1988 1940 A M DIMITS L L LODESTRO and D H E DUBIN Phys Fluids B 4 1992 274 T TAJIMA and F W PERKINS private communication 1983 T TAJIMA Computational Plasma Physics Addison Wesley Redwood City California 1989 M KOTSCHENRUETHER in Bull Am Phys Soc volume 33 p 2107 1988 S BRIGUGLIO G BETELLO G FOGACCIA G GRAZIADEI S C Guo C KAR F Ro MANELLI G VLAD and F ZONCA in Plasma Physics and Controlled Nuclear Fusion Research 1992 W rzburg volume II CN 56 D 2 3 pp 87 95 Vienna 1993 International Atomic Energy Agency S E PARKER and W W LEE Phys Fluids B 5 1993 77 76
71. ted kinetic energy of the I Fourier component K 3 IDIAGN 2 real and imaginary part of the myn at specific diagnostic radii as given in namelist diapos now IDIAGN 12 2 NRCHNL see Sect 3 2 Table 9 Quantities in the file ENERGY 32 This file is read by the plotting program epe3ak31 u and give global time evolution of the simulation from the point of view of the MHD module see Fig 4 An example of the input file xepe3ak31_datain for the program epe3ak31 u is listed hereafter ENERGY oo PRONF FPRORPF HFPOWHW o e 30 1 e 4 e 30 1 e 4 e 30 1 e 4 CHANGE TSTART IF 1 ADD IN THE NEXT LINE NEW TSTART CHANGE TEND IF 1 ADD IN THE NEXT LINE NEW TEND MODECH 1 MODES 2 ENERGY 3 PARAMETERS 4 MAGNETIC SIGNALS IC 1 MAGNETIC 2 KINETIC 3 TOTAL CHANGE MODE LIST IF 1 ENTER SEQUENCE OF MODE CHANGE CHANGE LIMITS IF 1 NEW LIMITS ARE ENTERED BY TERMINAL MODECH 1 MODES 2 ENERGY 3 PARAMETERS 4 MAGNETIC SIGNALS IC 1 MAGNETIC 2 KINETIC 3 TOTAL CHANGE MODE LIST IF 1 ENTER SEQUENCE OF MODE CHANGE CHANGE LIMITS IF 1 NEW LIMITS ARE ENTERED BY TERMINAL MODECH 1 MODES 2 ENERGY 3 PARAMETERS 4 MAGNETIC SIGNALS IC 1 MAGNETIC 2 KINETIC 3 TOTAL CHANGE MODE LIST IF 1 ENTER SEQUENCE OF MODE CHANGE CHANGE LIMITS IF 1 NEW LIMITS ARE ENTERED BY TERMINAL exit Note that in the above list only the sequence of input data suited for MODECH 1 have been shown dif
72. the input file xplot field input for the program plot field f is listed hereafter ENERGY PHIWRITE APWRITE TESTWRIK Te_vs_erg_DIII_D_i_interp txt 0 ips 0 no PostScript file 1 PS file 2 EPS name follows 30 char pippo eps 1 321 10 ifirst_step itot increm output time steps 1 1 phi 2 psi 0 phiplot toroidal angle for r theta plot only 1 iplO first plot 321 ipli last plot 0 0 l min l max Fourier components used 0 0 all 1 iflag rtheta 1 plot in r theta phiplot plane ph ap xxxx gif 0 iflag test superimposes to r theta phtest plot the i th test particle i iflag test 0 iflag_trajectory plot particle trajectory in rtest thtest phtest utest space trajec_ i _xxxx gif trajRZ_ i _xxxx gif trajXY_ i _xxxx gif 0 iflag_fourier_comp 1 plot fourier component profiles 0 iflag_contour 1 2D plot of Fourier components mn_ph ap_xxxx gif 2 contour plot mn_C_ph ap_xxxx gif 0 do both 0 iflag_power_spectrum 1 power power spectrum of field in the plane omega r 576 time window for Fourier transform 0 0 rO ri min max in r if 0 0 use max available interval 0 0 25 w0 wi min max in omega if 0 0 use max available interval 1 ihann Hanning function 0 off 1 on 3 ibuffer 1 no buffer n gt 1 zero buffer n 1 times O 0 001 ilog fac_zmin color scale 0 linear 1 log min value plotted 505 505 ndivx ndivy for power spectrum plot axes true logic_fill false only contour true fill
73. the right of the plot describing the m n poloidal and toroidal mode numbers of the curves are written every two modes because of space limitation although in the plot all the components are shown 36 4 3 Output file CTBO The CTBO file contains radial profile data of MHD quantities at the end of the simulation or each NCYCLE NSUBCY times if fortran is modified accordingly WRITE CHCTBO 1900 ETA VISC DT TIME REAL NR REAL LM REAL IRMODE WRITE CHCTBO 1900 RM L RN L L 21 LM WRITE CHCTBO 1900 R I I 1 NR 1900 FORMAT 6 E20 13 Cus WRITE CHCTB0 1950 PSI I L I 1 NP L 1 LM WRITE CHCTBO 1950 PHI I L I 1 NR L 1 LM CGV WRITE CHCTBO 1950 CUR I L I 1 NR L 1 LM WRITE CHCTBO 1950 W I L I 1 NR L 1 LM CGVKIN WRITE CHCTBO 1950 PREK I L I 1 NR L 1 LM CGVKIN CGV 1950 FORMAT 3 E20 13 E20 13 The variable IRMODE is an obsolete parameter included only for compatibility with old outputs This output file is read by the plotting program profilk u and gives the radial profile for each Fourier poloidal component of the various MHD variables the poloidal magnetic flux function PSI Wmn r the scalar potential PHI m n r the toroidal component of the current CUR jmn r A mn the toroidal component of the vorticity W Wm n r VAD ns and term proportional to the divergence of the energetic particle stress tensor which enters in the vor
74. ticity equation PREKm n r computed at the time t TIME see Fig 6 An example of the input file xprofilk_datain for the program profilk u is listed hereafter 0 0 title CTBO 1 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT 1 L mode index if 1 0 it plots from abs 1 to lm do not enter min max OEE X MIN X MAX Y MIN Y MAX if commas it takes computed values 1 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT 2 L M N if 1 lt 0 it plots from abs 1 to 1m do not enter min max 0 7133 X MIN X MAX Y MIN Y MAX if commas it takes computed values 8 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT 1 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT 3 L M N if 1 0 it plots from abs 1 to lm do not enter min max 8 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT 37 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT L M N if 1 lt 0 it plots from abs 1 to 1m do not enter min max 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT E 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT L M N if 1 lt 0 it plots from abs 1 to 1m do not enter min max 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT L M N if 1 lt 0 it plots from abs 1 to 1m do not enter min max 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT 1 PSI 2 W 3 PHI 4 J 5 Q 6 PREK 7 NEW TIME 8 BLANK PLOT L M N if 1 lt 0 it plots from abs 1 t
75. uced by assigning iflag_test 1 Right trajecto ries of the first test particle test test tests Utest produced by assigning iflag test 1 and flag trajectorysi 464a dh meu ela I2 16 42 10 11 12 13 Test particle trajectory projected on the poloidal cross section y 0 the cross indicates r 0 and on the equatorial plane Z 0 top view of the torus dotted line refers to r 0 the cross indicates the axis of symmetry of the torus produced by assigning iflag_test 1 and iflag EIA MIA EE ORE et de darn be See ef dn Left mn r produced by assigning iflag fourier comp 1 and iflag contour 0 Centre mn r contour plot and superimposed the curve m nq r Right frequency spectra of in the plane r w contour plot with superimposed the upper and lower Alfv n continua of the toroidal gap for a particular toroidal mode number automatically chosen from the first perturbed mode in this case n 2 produced by assigning itot 961 ipl0 ip11 481 iflag power spectrum 1 and other data from file xplot field input 4 5249 dai en Frequency spectra during the linear growth phase twao 144 and different values for twindow first row twindow 48 Wmin Wao 0 131 second row twindow 96 Wmin wao 0 0654 third row twindow 144 Wmin wao 0 0436 forth row twindow 192 wmin wao 0 0327 ihann and ibuffer first column ihann 0 ibuffer 1 second col umn ihann 1 ibuffer 1 third column ihann 1
76. ur namelists params diapos equipa paramk Then it contains a time sequence of some global fluid quantities namely WRITE CHENER 00003 TIMBUF JTBUF LM EZBUF JTBUF QOBUF JTBUF QABUF JTBUF 2 QABUF JTBUF 2 ENMODE JTBUF 1 1 ENBUF JTBUF I I 1 4 RM L RN CL ENMODE JTBUF K L K 1 IDIAGN 2 L 1 LM 00003 FORMAT F16 6 16 E16 6 3F10 6 amp 4E24 15 2F5 0 4E24 15 10X 4E24 15 amp 10X 4E24 15 10X 2E24 15 amp amp amp amp where the meaning of the quantities are listed in Table 9 The MHD module produces an output on the file ENERGY every NPRI time steps thus the output time step is given by Atoutput MHD DT NPRI in Alfv n time units 3l Quantities Comments TIMBUF simulation time t LM Fourier components considered in the simulation EZBUF electric field toroidal corrections neglected at r 1 significative only for nonlinear MHD simulations QOBUF q r 0 QABUF q r 1 2 QABUF internal inductance l toroidal corrections neglected ENBUF I I 1 volume integrated total magnetic energy I 2 volume integrated total kinetic energy I 3 resistive dissipation obsolete not corrected for toroidal terms I 4 viscous dissipation obsolete not corrected for toroidal terms RM L m 1 RN L n l ENMODE K L K 1 volume integrated magnetic energy of the l Fourier compo nent K 2 volume integra
77. utput time steps which will be read and available for plotting respectively frames will be produced every increm steps An output time step corresponds to NPRI simulation time steps i e to a sim ulation time Atoutput MHD DT NPRI in Alfv n time units time max is the limit of the abscissa in Alfv n time units if time max 0 the maximum of the abscissa will be calculated automatically by the input parameters of the run Note that for a simula tion with DT 0 02 NCYCLE 160 NSUBCY 3 NPRI 30 NOUT 20 we get a total number of time steps equal to NCYCLE NSUBCY NOUT 1 9601 t 0 is also counted corresponding to time max DT NCYCLE NSUBCY NOUT 192 0 and NCYCLE NSUBCY NOUT NPRI 1 321 output times y en min y en max are the limits of the ordinate if y en min 0 and or y en max 0 these limits are computed from the data iplO and ipli are the first and last time steps at which graphical frames will be produced in the above example only the frame number 161 will be produced corresponding to the normalized time twao 161 1 DT NPRI 96 The parameter 5 is a parameter required by the plotting routines HIGZ from CERN which identify the graphic window In Fig 5 an example of this plot is shown 35 Wiss Integr energy twy_ 96 00 it 3219 529 10 725 1 9 2 11 2 1 13 2 15 25 1 17 22 19 20 1 21 2 8 10 S 6 1 T Figure 5 Example of integrated total energy Wiot m n vs time The labels on
78. x ASPECT is the aspect ratio Ro a NERRE NR 1 RLOW is the radial coordinate of the mesh used by the MHD module PSIEQ are the radial profiles of the LMEQ equilibrium Fourier components of w MMODE 1 m l ANTOT is an obsolete quantity This file is used to produce movies of density Gy and ay local drive energetic particle profiles evolution see Fig 7 using the plot program plot density f An example of the input file xplot density input for the program plot density f is listed hereafter 40 ENERGY APWRITE OUTDNT 0 ips 0 no PostScript file 1 PS file 2 EPS name follows 30 char pippo eps 3 10 ifirst_step itot increm output time steps 0 1 plot n H density 2 beta_H 3 alpha_H 0 all 161 ipl0 first plot 161 ipli last plot 0 0 xldmin xldmax x axis if 0 0 use automatic values O 1 20 aldmin aldmax density if 0 0 use automatic values 0 0 011 albmin albmax beta H if 0 0 use automatic values 1 5 1 5 alamin alamax alpha H if 0 0 use automatic values 5 In the above example which refers to a simulation with DT 0 02 NCYCLE 160 NSUBCY 3 NWRITE 30 NOUT 20 we get a total number of time steps equal to NCYCLE NSUBCY NOUT 1 9601 t 0 is also counted corresponding to time_max DT NCYCLE NSUBCY NOUT 192 0 and NCYCLE NSUBCY NOUT NWRITE 1 321 output times Only the plots corresponding to the output time step iplO ipli 161 will be produced that is at the normalized time twao 16

HMGC User guide

Contents

Download Pdf Manuals

Related Search

Related Contents