Bers, A.

Fuchs, V.

Shoucri, M.

DSpace@MIT

PFC/JA-87-24THE APPLICATION OF THE FINITE ELEMENTLIBRARY TWODEPEP FOR THE NUMERICAL SOLUTIONOF THE RELATIVISTIC FOKKER-PLANCK EQUATIONbyM. Shoucri* V. Fuchs* A. BersJune 1987Plasma Fusion CenterMassachusetts Institute of TechnologyCambridge, Mass. 02139 U.S.A.*IREQ, Project Tokamak, C.P. 1000, Varennes, Quebec, CanadaJOL 2POSubmitted for publication in: Computer Physics CommunicationsTHE APPLICATION OF THE FINITE ELEMENT LIBRARYTWODEPEP FOR THE NUMERICAL SOLUTION OFTHE RELATIVISTIC FOKKER-PLANCK EQUATIONM. SHOUCRI, V. FUCHS, A. BERS*IREQ, Projet TokamakC.P. 1000, Varennes, Qufbec, Canada JOL 2POPermanent address:* Plasma Fusion Center, Massachusetts Institute of TechnologyCambridge, Mass. 02139, U.S.A.AbstractWe report on the application of a commercially available differentialequation solver library called TWODEPEP for the numerical solution of therelativistic Fokker-Planck equation. Results obtained by solving therunaway problem, and on the solution of the R.F. heating problem arepresented.2The interaction of radio-frequency waves with a plasma is described by aFokker-Planck equation with an added quasilinear diffusion term, and theresponse of a plasma to a homogeneous ohmic electric field is described bythe same Fokker-Planck equation with an advection term. There has been anextensive literature on numerical techniques to solve this equation whichis of high importance in the study of lower-hybrid RF-driven currents inTokamaks (see, for instance, the review article in Ref. [1]). In thisnote, we report on the application of a commercially available differentialequation solver library called TWODEPEP [2] for the numerical solution ofthe Fokker-Planck equation.This differential equation solver has been successfully applied in solvinga wide variety of problems of fluid dynamics, diffusion, heat transfer,etc. We have previously reported results obtained using thislibrary forthe numerical solution to the non-relativistic Fokker-Planck equation [3,5]. In this note we report on the application of this differential equa-tion solver for the numerical solution of the Fortran code used to solvethe relativistic Fokker-Planck equation. It illustrates how convenient isthe library for the problem we are studying, if we notice that the prepro-cessor Fortran code has no more than a few dozen lines. We note, however,that a physical quantities of interest plus graphic output, are obtainedfrom a separate code, which uses TWODEPEP results as input.In the model Fokker-Planck equation we consider the plasma is assumed homo-geneous, and azimuthally symmetric about the magnetic field. Relativisticeffects are included. The model is described in sufficient detail in3Refs. [1], [6], its characteristic feature being a linear collision opera-tor describing the interaction of fast test electrons with fixed Maxwellianfield electrons and ions. The electron distribution fonction f depends ontime and two momentum variables. In flux form the Fokker-Planck equationisbf ++tI V .S =0()where f is the test-electron distribution function and S the test electronflux in the presence of an electric field E and a quasi-linear diffusion*coefficient D:S = e f - D . grad f + Scoll (2)We use the normalized variablest + tv v + v/vt PP/Pt (3)where the thermal velocity unit vt is defined by the following relations[6]:vt2 1 V2 f(P)d3P = 5(1 - + 2(t = T- 2 e8(e e4where fm(P) denotes the relativistic Maxwellian distribution and 0 =Te/(me c 2). We will restrict ourselves to the case where 0 << 1 (forTe = 5 kV we have 0 = 0.01) so that the thermal velocity vt =(Te/me)/2. The thermal momentum is given byPt = mv tand the thermal collision frequency is given byv = 4un e 4lnA/m2vt3withvt (Te/m)1/2(5)(6)(7)The electric field is normalized to Eo = Ptv/e and the diffusioncoefficient is normalized to Do = vPt 2. Also, for 0 <(1, therelativistic Maxwellian distribution is given in our normalized units by:C..P2 /(Y+l)e (8)f (P) = 132m (2%)whereY = (I + P2p2 1/2Y =(l+Pth) (9)and5(10)Pth = vt/CIn the spherically symmetric coordinates (v, p), where p = Ps/P = cose,we have:+V.Sp(P2(P1 _ 6 (V 1 2 S )P oR(11)S= epf - Dpi (L a + (-S = (1- 2)i [ - ef + DP ) ! ) - (A(P) bf + F(P)) f(p + al)Here Z is the ion charge and the coefficients A(P), B(P) and F(P) are givenin our normalized units by:3A(P) = --T3 PB(P) = 213 2 PP , f (P ')0 2 d PP'Y'f m(P') dP')fP 3 p P4P2f (P')dP' - j -0 m 2P 3 Y (P') dP' +f P'y'f (P') dP')0where(12)(13)(14)(15)+ I(B(P) + Z)a62 PF(P) = f P'f (P') (3v' - v ' th dP' +3 P 2 m hdp+2P P'f t dP' (16)Y m tThe finite element code TWODEPEP [2] solves the 2-D partial differentialequationU Q QP5t = ax- + o-E+ R(17)Comparing Eq. (17) with Eq. (11) we can identify P + x, pt + y, U + f andQ = P2s 18)Q y= - P /(1-p12) (9The momentum-space domain in which the equation is solved is 0< P ' Pmaxand -1 <L I. The condition S = 0 at pL = ± 1 follows from Eq. (13),and represents an axisymmetric situation. At P = Pmax we always setf = 0, although one should impose a runaway flux condition Sp = Epf.With f = 0, however, the code will maintain required accuracy with consid-erably less computational effort in terms of required elements and memory.The f = 0 boundary condition at P = Pmax will distort the solution onlyin a narrow boundary layer close to Pmax. To avoid the singularity at7P = 0, we impose the boundary condition f = fm (fm is defined in Eq.(8)) at P = Pmin (Pmin could be made very small, say 0.01). The impo-sition of f at P = Pmin will slightly violate density conservation in thebulk when e and/or D are non-zero, but this effect is very small since forthe values of e and D considered the number of particles transferred to thetail is orders of magnitude smaller than the number of particles in thebulk. Thus when e * 0, an imposed fixed boundary condition in the bulkrepresents an equilibrium situation in which electrons are being fed intothe distribution at the same rate at which they run away. Hence we assumesteady-state and solve for the equilibrium:bQ bQb-2E+ --l = 0 (20)x yThe functions Qx and Qy are defined in Eqs. (18) and (19) (see subrou-tines FXX and FXY in the Appendix, in which U + f, UX + and UY + ).tion (20) is a linear equation for f and TWODEPEP requires only one itera-tion to calculate the solution (parameter TF = 1 in the main, and NOUT = 1denotes that the output should be given after one iteration). Also in themain NEQ = 1 indicates we are solving one equation, NTF = 9300 denotes thenumber of elements we are using and NDIM = 2 is a memory storage param-eter. With these parameters the code would require a memory close to 5.5million bytes on an IBM machine using double precision arithmetic. OXX andOXY are the functions Qx and Qy (defined in Eqs. (18) and (19)) whichare calculated in the functions FXX and FXY respectively. D3EST is a func-tion which distributes the elements in such a way as to have more elementsin regions where more precision is required. UPRINT is the output solution8given by U ~ f. The output solution is given on a grid defined by x = XA +i * Hx ( i = 0, ... NX) and y = YA + j * HY (j = 0, ... NY). MWR = 8 isthe file number on which the results are stored, and the INTEGRAL functionintegrates a function of the solution over the domain where the solution iscalculated. In the Appendix presented, we are calculating the rf powerabsorbed defined by:1 oP = -2% J d f D (P3/y) dP (21)d aoThe remaining parameters in the main define the boundary conditions. Formore details see Ref. [2]. Finally, the calculations of the integrals inEqs. (14) - (16) has been simplified by noting that for small values of P,and for the small values of H we are considering, the quantity v = P/Y = vand the integrals in Eqs. (14) - (16) can be calculated from well-knownerror functions integrals [1], while for large P, we have the asymptoticexpressions:1A(P) =- (22)vB(P) = - (1 - -2-) (23)2v vF(P) = -1 (24)vWe present in the following two sample results obtained using the programin the Appendix.9a) The runaway distribution function:The results presented in this case are with e = 0.04, D = 0 and Pth= 0.0442(i.e. 1 kV plasmas). In Fig. (1), we plot the perpendicular velocitymoment of the distribution function f, defined byMF = 2% f dPOf (25)0as a function of PH. The function F is the parallel distribution func-tion. It shows clearly the formation of a runaway tail. We show inFig. (2) perpendicular cuts of f, indicating the presence of a hot tailpopulation, a result similar to what has been recently presented for thenon-relativistric case [5]. The flux FR of electrons through the surfaceof a sphere of radius R is defined by3f=f 3 +FR fdPf= d PdivS = 2nP2 f dpS, (26)For the present case (e = 0.04) the result for FR is 0.14 x 10- . Thisis in very good agreement with the value 0.13 x 10-5 obtained from theexpression [7]:2 1 2 1/2rR =rNR ex th 2 + -=/2 (l + Z) )}(27)10where rNR is the runaway production rate for the non-relativistic case(rNR = 0.196 x 105 for e = 0.04). Table I summarizes some of theresults obtained for different values of e for a 1 kV plasma (@th= 0.0442).JR denotes the relativistic current (in units of envt)JR vI fd3P (28)calculated for different values of the normalized electric field e and Z.rNR and rR are the runaway production rates calculated for the non-relativistic case and the relativistic case respectively and in the lastcolumn TR is calculated from Eq. (27), showing a very good agreement withthe value in the previous column. Finally, we show in Fig. (3) contour-plots for the 2-D solution f obtained for e - 0.04.b) Lower-Hybrid Current Drive.We consider the case when e = 0 and the quasilinear diffusion term is givenby:3 for vj <P /Y < v2 (29)0 otherwiseD * 0 inside a region bounded by the hyperbolae whose equations are givenbyP2 12 P111 2__2 9th th = 1 (30)vi, 211Figure (4) shows a contour plot of f obtained for the case vi = 4, v2 = 7and @th = 0.1. The region where D * 0 is clearly shown to be bounded bythe hyperbola whose equations is given by Eq. (30). Figure (5) shows aplot of the parallel distribution function F as defined in Eq. (25). Aplateau formation is apparent in the region where D * 0. The currentassociated with this distribution is JR = 0.45 x 10-2 and the powerabsorbed as defined in Eq. (21) is 1.57 x 10-4 resulting in a figure ofmerit JR/Pd = 28.To conclude, using TWODEPEP library and a Fortran code of no more than afew dozen lines, we have obtained results for the numerical solution of thelinearized Fokker-Planck equation which compares favourably with what hasbeen previously reported in the literature for the numerical solution ofthis problem [1, 6]. We are now in the process of extending these resultsto higher plasma temperature and to velocity spectra extending to higherrelativistic speeds. In this latter case, however, we have been limited bythe computer memory which for the present runs is close to 5.5 millionbytes. Typical runs would take 30 min CPU time on the IBM 3084. Higherrelativistic speeds will require a larger number of elements and conse-quently also a larger memory for the library. The use of cubic elementswill require an even larger memory.12TABLE Ie z 102 NR R R (EQ. 27)0.005 1 1.9 0.76 x 10-30 0.46 x 10-38 0.3 x 10-360.01 1 3.8 0.13 x 10-16 0.18 x 10-18 0.18 x 10-0.03 1 12.7 0.83 x 10-7 0.5 x 10-7 0.45 x 10-70.04 1 18.8 0.196 x 10-5 0.14 x 10. 5  0.13 x 10-50.04 4 8.25 0.158 x 10-6 1.02 x 10-7 0.94 x 10-713REFERENCES1. C.F.F. Karney, "Fokker-Planck and Quasilinear Codes," Report PPPL-2290(Princeton), 1985.2. TWODEPEP, IMSL Problem-Solving Software System for PartialDifferential Equations (IMSL, Houston 1983), IMSL No. TDP-0005,Edition 5.3. V. Fuchs, R.A. Cairns, M.M. Shoucri, K. Hizanidis and A. Bers, Phys.Fluids 28, (1985) 3619.4. V. Fuchs, R.A. Cairns, C.N. Lashmore-Davies and M.M. Shoucri, Phys.Fluids 29, (1986) 2931.5. M. Shoucri, V. Fuchs and R.P. Gupta "Numerical Solution of a Two-Dimensional Fokker-Planck Equation" Report IREQ, Varennes (Nov. 1983)M. Shoucri, V. Fuchs, K. Hizanidis and A. Bers "Enhancement ofElectron Runaway and Toroidal Current in the presence of Lower HybridWaves" Report Tokamak de Varennes, TRI 165 (December 1984); V. Fuchs,M. Shoucri, A. Bers and R.A. Cairns "Lower Hybrid Current Drive in thePresence of an Ohmic Electric Field" Report Tokamak de Varennes TVRI187 (July 1985).6. C.F.F. Karney and N.J. Fisch, Phys. Fluids 28, (1985) 116.7. J.W. Connor, R.J. Hastie, Nucl. Fusion 15, (1975) 415.14Appendix0XXOX Y03E3TUPRINTWXNYH XH YX AYANOUT1 9300 2FX.X (X r Y rU UX iUY 9 FCT)FXY(XFYPUUXPUY)WEIGHT(XtY)Usoo0.300.04Function Qx in Eq.(48)- Function Oy in Eq. (49)-1.1TF 1MWR aINTEGRALFNCT(XY'UUXUYNOUTFT)ARC-1FB1 FB(XPY)ARC -3FB1 0.ARC =ZGB1 0.ARC-4GB1 0.XGRID 0.1,60.1YGRID -1.,1.IX -IF-3IY 2,4ADD.Define boundarycurve andboundary conditionsDOUBLE PRECISION FUNCTION FXX(XYUUXFUYPFCT)IMPLICIT REALX8(A-HO-Z)COMMON/01/ZCPIEPSDOVP1,VPZFACTBTHDENTAUGAMA=DSQRT(1.+XXNBTHZBTH)VEL=X/GAMAS=VEL/DSORT(2.DO)VPAR=XxY/GAMAD=DOIF(VPAR.LT.VP1.OR.VPARGT.VP2) D=0.DOSD=D(YxUX+(1.-YxY)XUY/X)SE=EPSzUFU=(SxDEXP(-SxS)+ZxSxDEXP(-SESRCXC)/C)X2./DSQRT(3.1416D0)IF(X.LE.3.2)$BU=DERF(S)-FU+ZxDERF(CXS)/CxxZIF(X.GT.3.2) BU-i.A A-BU/VELxx3FF=BU/VELxZ2FXX=XxXxYx(SE-SD)-(AAxUX+FFXU)xXRX - Define Ox in Eq. (48)15FCT=FXXXX/GAMARETURNENDDOUBLE PRECISION FUNCTION FXY(XYPUPUXPUY)IMPLICIT REALx8(A-HO-Z)COMMON/D1/Z.CPIEPSDOVPIVP2,FACTBTHDENTAUGAMA=DSORT (1.+AX zXBTHxBTH)VEL=X/GAMAS=VEL/DSQRT(2.DO)VPAR=XxY/GAMAD=DOIF(VPAR.LTVPI.OR.VPAR.GT.VP2) D=O.DOSD=Dx(YxUX+(1.-YxY)ZUY/X)SE=EPS*UFU=(SxDEXP(-SIS)XFACT$ +ZRSXDEXP(-SMSXCC)/C)*2./DSRT(3141600)FU2=14-1./(2.XSXS)FU2C=l.-./(2,XSXSZCXC)IF(X.LE.3.2)$AU=DERF(S)*FU2xFACT+ZROERF(CXS)XFUZC+FU/(2.xSXS)IF(X.GT.3.2) AU=1.-1./VELxxZ2+ZBB=AU/(Z.xVEL)FXY -Xx(1.-YxY)(BBxUY/X-SE+SD) - Define Qy in Eq. (49)RETURNENDDOUBLE PRECISION FUNCTION FNCT(XYUUXPUYN~UTPlT)IMPLICIT REALXS(A-HO-Z)COMMON/Dl/ZCPI.EPSD0,VP1,VPZFACTBTHDENTAUFNCT=O.DONOU=NOUT+1NOUM=NOU+1GAMA=DSQRT(I.+XxXxBTHxBTH)VPAR=XY/GAMA Function toVELX/TAMA be integratedFU=(SXDEXP(-SxS)+ZESxDEXP(-S*SxCXC)/C)x2./DSQRT(3.141600) at the endIF(X.LE.3.2) BU=DERF(S)-FU+ZXDERF(CxS)/Cx*2IF(X.GT.3.2) BU=1.AA=BU/VELxx3FF=BU/VELxx2D=DOIF('PARLT,VP1.OR.VPARGT.VP2) D=0.D0FNCT=-DzXxYx(YxUX+(1,-YxY)xUY/X)x2.x3.14l6xXx/GAMARETURN16ENDDOUBLE PRECISION FUNCTION FB(XPY)IMPLICIT REALX8(A-HPO-Z) Define boundaryCOMMON/D1/ZCPIEPS'DOVPIVP2,FACTPBTHYDENTAU condition forGAMA=DSQRT(i. +XxX*BTHxBTH) ARC= -4FB=0,0634936xDEXP (-XxX/(1,+GAMA))RETURNENDDOUBLE PRECISION FUNCTION WEIGHT(XPY)IMPLICIT REALXS.(A-HO-Z)COMMON/Di/ZCPIEPS,DOVPIVP2'FACTBTHDENTAUVi=VPi-i.V2=VP2+1#IF(EPS.NE.0.) VCRIT=DSQRT(1./EPS)GAMA=DSORT(1.+XXXZBTHXBTH)VPAR=XxY/GAMAWDO=2.WE0=Z.W0-0 .WE=0.IF(EPS.NE.0..AND.VPAR.GT.VCRIT) WE=WEOIF(DO.NE.0..AND.(VPAR.GT.Vl.AND.VPAR.LT.VZ)) WD=WDOWEIGHT=i.+WD+WERETURNENDBLOCK DATAIMPLICIT REALX'(A-HPO-Z)COMMON/D1/ZCpPIEPSDOVPiVP2,FACTBTHDEN.TAUDATA Z/1.DO/,C/42.900/,PI/3.14159265D0/,FACT/1.ODO/DATA EPS/0.00/,00/3.ODO/,VPI/4.0000/,VP2/16.0D0/,BTH/0.054DO/DATA DEN/0.5D0/,TAU/15.0D-3/ENDEND.17FIGURE CAPTIONSFigure 1 The perpendicular velocity moment of the distribution function f,defined by Eq. (25), for the runaway problem.Figure 2 Perpendicular cuts of the distribution function f at differentvalues of P11, for the runaway problem.Figure 3 Contourplots of the distribution function f, for the runawayproblem.Figure 4 Contourplots of the distribution function f for the RF problem.Figure 5 The perpendicular velocity moment of the distribution function f,defined by Eq. (25), for the RF problem.91-30 -20 -10 0pit40 20 3018CC.94-C)C'4-C04-1~4-Li,0I II I4236Distribution functiono o on0 0oc000-000000"I000000I-0000U'0000r.30-00I-r.3',.I00000019-DDistribution function- ~b 0 *0 b 0bItDistribution funtiono 0 0 0 0 0 0-O--rOIt 0 OAcJ cJ200090--0--.0U)j %,C>2104V3018A-30 -20 -40 0P1i10 20 3022C0C-)CC0~01~(j~05-209n R

The application of the finite element library twodepep for the numerical solution of the relativistic fokker-planck equation

http://dspace.mit.edu/bitstream/1721.1/94951/1/87ja024_full.pdf

The application of the finite element library twodepep for the numerical solution of the relativistic fokker-planck equation

Abstract

Similar works

Full text

Available Versions

DSpace@MIT