1- # ---
2- # title: Tokamak profile
3- # id: demo_tokamak_profile
4- # date: 2023-04-20
5- # author: "[Tiancheng Liu](https://github.com/TCLiuu)"
6- # julia: 1.9.0
7- # description: Tracing passing and trapped proton in a Tokamak
8- # ---
9-
10- # This example shows how to trace protons in a stationary magnetic field that
11- # corresponds to an ITER-like Tokamak.
12-
13- import DisplayAs # hide
14- using TestParticle
15- using TestParticle: getB_tokamak_profile
16- using OrdinaryDiffEq
17- using StaticArrays
18- using GeometryBasics
19- using CairoMakie
20- CairoMakie. activate! (type = " png"
21-
22- # # parameters from ITER, see http://fusionwiki.ciemat.es/wiki/ITER
23- const R₀ = 6.2 # Major radius [m]
24- const Bζ0 = 5.3 # toroidal field on axis [T]
25- const a = 2.0 # Minor radius [m]
26-
27- # # variable must be a radius normalized by minor radius.
28- function q_profile (nr:: Float64 )
29- return nr^ 2 + 2 * nr + 0.5
30- end
31-
32- function B (xu)
33- SVector {3} (getB_tokamak_profile (xu[1 ], xu[2 ], xu[3 ], q_profile, a, R₀, Bζ0))
34- end
35-
36- function E (xu)
37- SA[0.0 , 0.0 , 0.0 ]
38- end
39-
40- # Passing proton in a Tokamak
41-
42- # # initial velocity for passing particle
43- v₀ = [0.0 , 2.15 , 3.1 ] .* 1e6
44- # # initial position, [m]. where q≈2, (2, 1) flux surface.
45- r₀ = [7.3622 , 0.0 , 0.0 ]
46- stateinit = [r₀... , v₀... ]
47-
48- param = prepare (E, B; species= Proton)
49- tspan = (0.0 , 4e-5 ) # [s]
50-
51- prob = ODEProblem (trace!, stateinit, tspan, param)
52- sol = solve (prob, Vern7 (); dt= 2e-11 )
53-
54- fig1 = orbit (sol)
55- # # contruct the topology of Tokamak
56- nθ = LinRange (0 , 2 π, 30 )
57- nζ = LinRange (0 , 2 π, 30 )
58- nx = [R₀* cos (ζ) + a* cos (θ)* cos (ζ) for θ in nθ, ζ in nζ]
59- ny = [R₀* sin (ζ) + a* cos (θ)* sin (ζ) for θ in nθ, ζ in nζ]
60- nz = [a* sin (θ) for θ in nθ, ζ in nζ]
61- points = vec ([Point3f (xv, yv, zv) for (xv, yv, zv) in zip (nx, ny, nz)])
62- faces = decompose (QuadFace{GLIndex}, Tesselation (Rect (0 , 0 , 1 , 1 ), size (nz)))
63- tor_mesh = GeometryBasics. Mesh (points, faces)
64- # # plot the surface of Tokamak
65- wireframe! (fig1[1 , 1 ], tor_mesh, color= (:blue , 0.1 ), linewidth= 0.5 , transparency= true )
66- Label (fig1[1 , 1 , Top ()], " passing particle" = (0 , 0 , 10 , 0 ), fontsize = 22 )
67-
68- fig1 = DisplayAs. PNG (fig1) # hide
69-
70- # Trapped proton in a Tokamak that shows the banana orbit
71-
72- # # initial velocity for trapped particle
73- v₀ = [0.0 , 1.15 , 5.1 ] .* 1e6
74- # # initial position, [m]. where q≈1, (1, 1) flux surface.
75- r₀ = [6.6494 , 0.0 , 0.0 ]
76- stateinit = [r₀... , v₀... ]
77-
78- param = prepare (E, B; species= Proton)
79- tspan = (0.0 , 4e-5 )
80-
81- prob = ODEProblem (trace!, stateinit, tspan, param)
82- sol = solve (prob, Vern7 (); dt= 1e-11 )
83-
84- fig2 = orbit (sol)
85-
86- wireframe! (fig2[1 , 1 ], tor_mesh, color= (:blue , 0.1 ), linewidth= 0.5 , transparency= true )
87- Label (fig2[1 , 1 , Top ()], " trapped particle" = (0 , 0 , 10 , 0 ), fontsize = 22 )
88-
89- # # banana orbit
1+ # ---
2+ # title: Tokamak profile
3+ # id: demo_tokamak_profile
4+ # date: 2023-04-20
5+ # author: "[Tiancheng Liu](https://github.com/TCLiuu)"
6+ # julia: 1.9.0
7+ # description: Tracing passing and trapped proton in a Tokamak
8+ # ---
9+
10+ # This example shows how to trace protons in a stationary magnetic field that
11+ # corresponds to an ITER-like Tokamak.
12+
13+ import DisplayAs # hide
14+ using TestParticle
15+ using TestParticle: getB_tokamak_profile
16+ using OrdinaryDiffEq
17+ using StaticArrays
18+ using GeometryBasics
19+ using CairoMakie
20+ CairoMakie. activate! (type = " png"
21+
22+ # # parameters from ITER, see http://fusionwiki.ciemat.es/wiki/ITER
23+ const R₀ = 6.2 # Major radius [m]
24+ const Bζ0 = 5.3 # toroidal field on axis [T]
25+ const a = 2.0 # Minor radius [m]
26+
27+ # # variable must be a radius normalized by minor radius.
28+ function q_profile (nr:: Float64 )
29+ return nr^ 2 + 2 * nr + 0.5
30+ end
31+
32+ function B (xu)
33+ SVector {3} (getB_tokamak_profile (xu[1 ], xu[2 ], xu[3 ], q_profile, a, R₀, Bζ0))
34+ end
35+
36+ function E (xu)
37+ SA[0.0 , 0.0 , 0.0 ]
38+ end
39+
40+ # Passing proton in a Tokamak
41+
42+ # # initial velocity for passing particle
43+ v₀ = [0.0 , 2.15 , 3.1 ] .* 1e6
44+ # # initial position, [m]. where q≈2, (2, 1) flux surface.
45+ r₀ = [7.3622 , 0.0 , 0.0 ]
46+ stateinit = [r₀... , v₀... ]
47+
48+ param = prepare (E, B; species= Proton)
49+ tspan = (0.0 , 4e-5 ) # [s]
50+
51+ prob = ODEProblem (trace!, stateinit, tspan, param)
52+ sol = solve (prob, Vern7 (); dt= 2e-11 )
53+
54+ fig1 = orbit (sol)
55+ # # contruct the topology of Tokamak
56+ nθ = LinRange (0 , 2 π, 30 )
57+ nζ = LinRange (0 , 2 π, 30 )
58+ nx = [R₀* cos (ζ) + a* cos (θ)* cos (ζ) for θ in nθ, ζ in nζ]
59+ ny = [R₀* sin (ζ) + a* cos (θ)* sin (ζ) for θ in nθ, ζ in nζ]
60+ nz = [a* sin (θ) for θ in nθ, ζ in nζ]
61+ points = vec ([Point3f (xv, yv, zv) for (xv, yv, zv) in zip (nx, ny, nz)])
62+ faces = decompose (QuadFace{GLIndex}, Tesselation (Rect (0 , 0 , 1 , 1 ), size (nz)))
63+ tor_mesh = GeometryBasics. Mesh (points, faces)
64+ # # plot the surface of Tokamak
65+ wireframe! (fig1[1 , 1 ], tor_mesh, color= (:blue , 0.1 ), linewidth= 0.5 , transparency= true )
66+ Label (fig1[1 , 1 , Top ()], " passing particle" = (0 , 0 , 10 , 0 ), fontsize = 22 )
67+
68+ fig1 = DisplayAs. PNG (fig1) # hide
69+
70+ # Trapped proton in a Tokamak that shows the banana orbit
71+
72+ # # initial velocity for trapped particle
73+ v₀ = [0.0 , 1.15 , 5.1 ] .* 1e6
74+ # # initial position, [m]. where q≈1, (1, 1) flux surface.
75+ r₀ = [6.6494 , 0.0 , 0.0 ]
76+ stateinit = [r₀... , v₀... ]
77+
78+ param = prepare (E, B; species= Proton)
79+ tspan = (0.0 , 4e-5 )
80+
81+ prob = ODEProblem (trace!, stateinit, tspan, param)
82+ sol = solve (prob, Vern7 (); dt= 1e-11 )
83+
84+ fig2 = orbit (sol)
85+
86+ wireframe! (fig2[1 , 1 ], tor_mesh, color= (:blue , 0.1 ), linewidth= 0.5 , transparency= true )
87+ Label (fig2[1 , 1 , Top ()], " trapped particle" = (0 , 0 , 10 , 0 ), fontsize = 22 )
88+
89+ # # banana orbit
9090fig2 = DisplayAs. PNG (fig2) # hide
0 commit comments