ITER Equilibrium in Cylindrical Coordinates

Load Dependencies

using CairoMakie
using GeometricIntegrators
using LaTeXStrings
using LinearAlgebra

using ChargedParticleDynamics
using ChargedParticleDynamics: InitialConditions, charged_particle, guiding_center, pauli_particle
using ChargedParticleDynamics: md

import ChargedParticleDynamics.ChargedParticle3d
import ChargedParticleDynamics.GuidingCenter3d
import ChargedParticleDynamics.GuidingCenter4d
import ChargedParticleDynamics.PauliParticle3d

import .ChargedParticle3d.TokamakIterCylindrical: from_cartesian, g₁₁, g₂₂, g₃₃
import .ChargedParticle3d.TokamakIterCylindrical: R₀, B, b, aₚ, bₚ, cₚ, b⃗, ḡ, DF̄, J

Magnetic Field

fig, ax = plot_fieldlines(ChargedParticle3d.TokamakIterCylindrical; xrange = (3., 11.), yrange = (-4., +4.),  levels = 25)
fig

Initial conditions and solver options

Set and initialize initial conditions:

X₀ = from_cartesian(0, [7.0, 0, 0])
E₀ = 1E6
θ₀ = 0.
α₀ = π*5/16
m₀ = md

ics = InitialConditions(X₀, θ₀, α₀, E₀, m₀, 1, aₚ, bₚ, cₚ, b⃗, B, ḡ, DF̄, J)

Charged Particle Initial Conditions with
  x   = [7.0, -0.039356746343052945, 0.0005129847414674947]
  X   = [7.0, 0.0, 0.0]
  ρ   = [3.220094336611642e-19, -0.039356746343052945, 0.0005129847414674947]
  v   = [-0.1862898084737484, 0.007620581643085356, 0.01193182240348234]
  v∥  = [0.0, 0.007620581643085356, 0.01193182240348234]
  v⟂  = [-0.1862898084737484, 0.0, 0.0]
  |v| = 0.20429884217319663
  u   = 0.0838696856565162
  μ   = 0.003681107352162368
  θ   = 0.0
  α   = 0.9817477042468103
  ω   = 2.258748612191711e8
  M   = 3.3435837724e-27
  E   = 1.0e6
  C   = 1

Obtain charged particle initial conditions:

x₀, v₀ = charged_particle(ics)
p₀ = ChargedParticle3d.TokamakIterCylindrical.charged_particle_3d_pᵢ(0, x₀, v₀)

3-element Vector{Float64}:
 -0.2786657145332284
 -1.9863382989370268
  1.7868148787681433

Obtain Pauli particle initial conditions:

X₀, u₀, μ = pauli_particle(ics)

([7.0, 0.0, 0.0], [0.0, 0.007620581643085356, 0.01193182240348234], 0.003681107352162368)

Obtain guiding center initial conditions:

q₀, μ = guiding_center(ics)

([7.0, 0.0, 0.0, 0.0838696856565162], 0.003681107352162368)

The 3D model can use the guiding centre initial condition as it stands. It used to need displacing off the midplane by sqrt(eps()), because the only constraint pair it implemented divides by b₁ = b_R, which vanishes there; TokamakIterCylindrical now defaults to the pair that divides by b₃ = b_φ instead.

q3₀ = q₀

4-element Vector{Float64}:
 7.0
 0.0
 0.0
 0.0838696856565162

Parameters:

parameters = (μ=μ,)

(μ = 0.003681107352162368,)

Time step size and integration intervals:

Δt1 = 0.01
Δt2 = 0.1
Δt3 = 1.0
Δt4 = 10.0
Δt5 = 50.0
timespan = (0., 10_000.)
timespan_long = (0., 1_000_000.);

(0.0, 1.0e6)

Nonlinear solver options. f_abstol is an absolute bound on the residual and has to stay above its round-off floor, which for these ITER-scale models is ‖ϑ‖ eps ≈ 3.3E-15. Below that the Newton iteration has no reachable stopping criterion and runs to its limit on nearly every step. The library default does not clear the floor either: GeometricIntegratorsBase scales it as max(8, solversize(method, problem)) * eps(datatype(problem)), which is 1.8E-15 here. See Model Audit for the full reasoning:

options = (f_abstol = 1E-12, max_iterations = 50, warn_iterations = 50)

(f_abstol = 1.0e-12, max_iterations = 50, warn_iterations = 50)

Charged Particle

Integrate charged particle dynamics and convert result to cartesian coordinates:

code = ChargedParticle3d.TokamakIterCylindrical.podeproblem(x₀, p₀; timestep=Δt1, timespan=timespan)
csol = integrate(code, PartitionedGauss(1); options...)
ccar = cartesian_solution(csol, ChargedParticle3d.TokamakIterCylindrical);

(t = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09  …  9999.91, 9999.92, 9999.93, 9999.94, 9999.95, 9999.96, 9999.97, 9999.98, 9999.99, 10000.0], R = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.998138186966312, 6.996280612941102, 6.99443141068637, 6.992594697447767, 6.99077456577642, 6.988975074378906, 6.987200239016345, 6.985454023473781, 6.983740330621002  …  6.22059111990751, 6.221896787221355, 6.223281118534398, 6.224740349126757, 6.226270505335941, 6.22786741618185, 6.2295267255392295, 6.231243904817405, 6.233014266105845, 6.23483297574298], X = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[6.999999078963313, 6.998136785350421, 6.996278624065781, 6.9944287249548625, 6.992591202363984, 6.990770145973562, 6.988969611659651, 6.98719361240473, 6.985446109278808, 6.9837310025119965  …  -6.005972977616585, -6.007133259497455, -6.00837163635428, -6.009684380907654, -6.011067558392678, -6.0125170380738835, -6.01402850530365, -6.015597474084324, -6.01721930009295, -6.018889194126462], Y = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.003590893032779872, 0.004429153576905424, 0.005275363091232647, 0.006129463468924294, 0.006991380151815506, 0.007861022279658587, 0.00873828287348913, 0.009623039052801112, 0.010515152286145569, 0.011414468674694403  …  -1.619889648467311, -1.620601627012295, -1.6213260498000222, -1.622062284846451, -1.6228096666138516, -1.6235674979534576, -1.6243350521375373, -1.6251115749739653, -1.625896286996157, -1.6266883857210526], Z = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[-0.039356746343052945, -0.03932433269939358, -0.039379419304171714, -0.039521735187494085, -0.03975081572982512, -0.04006600340946461, -0.04046644898399032, -0.04095111310551217, -0.041518768368548016, -0.04216800178827693  …  -0.9076505916711205, -0.909173140769324, -0.9106267426367134, -0.9120073166013223, -0.9133109906801823, -0.9145341123734546, -0.9156732587876171, -0.9167252460601127, -0.9176871380604796, -0.9185562543456504])

Plot solution and energy error:

fig, ax = plot_trajectory_poloidal(ccar.R, ccar.Z, ChargedParticle3d.TokamakIterCylindrical; linewidth = 1)
fig

fig, ax = plot_trajectory_3d(ccar.X, ccar.Y, ccar.Z)
fig

cham = compute_invariant(csol.t, csol.q, csol.p, code.parameters, ChargedParticle3d.TokamakIterCylindrical.hamiltonian)
plot_invariant_error(csol.t, cham; label="H", padding_right=30)

Integrate charged particle dynamics with large time step and convert result to cartesian coordinates:

code_Δt2 = ChargedParticle3d.TokamakIterCylindrical.podeproblem(x₀, p₀; timestep=Δt2, timespan=timespan)
csol_Δt2 = integrate(code_Δt2, PartitionedGauss(1); options...)
ccar_Δt2 = cartesian_solution(csol_Δt2, ChargedParticle3d.TokamakIterCylindrical);

(t = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9  …  9999.1, 9999.2, 9999.3, 9999.4, 9999.5, 9999.6, 9999.7, 9999.8, 9999.9, 10000.0], R = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.982358403006357, 6.968458106569868, 6.961256929176343, 6.962300181903745, 6.971381026610796, 6.986586109319926, 7.004715117709091, 7.021971626154481, 7.034764887447843  …  7.279018818995169, 7.2790521507651205, 7.286841385637905, 7.30088878857537, 7.318467693158941, 7.336172801705416, 7.350593918866171, 7.358975878323361, 7.3597410479066525, 7.352784287761649], X = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[6.999999078963313, 6.982347558404731, 6.96842440009305, 6.961185838691547, 6.962177400674671, 6.971194473459793, 6.986327374417188, 7.0043797480727275, 7.02155810887263, 7.034272666598382  …  6.885182044795722, 6.883031242971213, 6.888202489911516, 6.89913294970817, 6.913128094458155, 6.9269056406931115, 6.9372450612015975, 6.941610696815994, 6.938631951071888, 6.928358053715153], Y = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.003590893032779872, 0.012306164023587958, 0.02167402148006746, 0.03146034338516039, 0.041348089094939755, 0.051000307770133826, 0.06012720192379662, 0.06854359963501414, 0.07620526388401318, 0.08321702719965376  …  -2.3618601095976097, -2.3682231993290848, -2.3771251202707275, -2.3885019667977216, -2.4017971033008845, -2.4160731821317505, -2.4301979180524302, -2.443065064953214, -2.453808252452998, -2.461969021306922], Z = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[-0.039356746343052945, -0.04274895253676522, -0.05358348898469442, -0.06941987350642155, -0.08674601808704696, -0.10172914167767044, -0.11104053067920731, -0.11256456482690617, -0.1058348486836121, -0.09211236563995544  …  0.23338119887026826, 0.2145379872932519, 0.1973853546445848, 0.18508379741896175, 0.1798255875891623, 0.18240984351384537, 0.19210138465564833, 0.20678967683170363, 0.2234057380466143, 0.2385126738306136])

Pauli Particle with Symplectic Integrator

Integrate charged particle dynamics and convert result to cartesian coordinates:

hode = PauliParticle3d.TokamakIterCylindrical.hodeproblem(X₀, u₀; parameters = parameters, timestep=Δt3, timespan=timespan)
hsol = integrate(hode, PartitionedGauss(1); options...)
hcar = cartesian_solution(hsol, PauliParticle3d.TokamakIterCylindrical);

(t = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …  9991.0, 9992.0, 9993.0, 9994.0, 9995.0, 9996.0, 9997.0, 9998.0, 9999.0, 10000.0], R = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 7.000218761141969, 6.9999999117627345, 6.999682716546662, 6.999671599626917, 6.999066237995943, 6.998711774229864, 6.998335908170604, 6.997471736575934, 6.9970398765712165  …  5.83847503837543, 5.839237543730642, 5.840100721065733, 5.8406440022063375, 5.841713597229706, 5.84241331593058, 5.843291051652109, 5.844390660043766, 5.845095810206043, 5.846321252513696], X = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.9997211061009805, 6.99801031707669, 6.99520559645081, 6.9917137816740125, 6.986635998983829, 6.980814527465681, 6.9739846838298964, 6.965677732307828, 6.956813247226993  …  -5.641146036642602, -5.643780389579486, -5.646598521258317, -5.649194556111499, -5.652367495636069, -5.655274368669487, -5.658421990399065, -5.66185968035947, -5.66499716509576, -5.668697673732491], Y = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.08346939947350995, 0.1668842915510738, 0.2503133949778289, 0.3336778353042345, 0.41694798423535007, 0.5001949938452352, 0.5833036201636984, 0.6662915531581417, 0.7492105695562288  …  1.5050787909622656, 1.4981448549160683, 1.490872755105504, 1.4831463945642962, 1.4752489438516025, 1.4668555379562394, 1.458187536939572, 1.4492229459702088, 1.4397750344915843, 1.4301534433105962], Z = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.00822547311397348, 0.016820456442619314, 0.024899349641928747, 0.03331836604044591, 0.04177155110327058, 0.04982403243743794, 0.05838201203210591, 0.06661932922032276, 0.07477324893333392  …  -0.8936857913199369, -0.8935995955901991, -0.8930377620831156, -0.8928316577731679, -0.8925975431840649, -0.8921006127625265, -0.8920896088981086, -0.8916341043419427, -0.8914321801092542, -0.8913281418473084])

Plot solution and energy error:

fig, ax = plot_trajectory_poloidal(ccar.R, ccar.Z, ChargedParticle3d.TokamakIterCylindrical; label="Charged Particle Δt=0.01", linewidth = 1)
plot_trajectory_scatter!(fig, ax, hcar.R, hcar.Z; label="Pauli Particle Δt=1.0")
axislegend(ax)
fig

fig, ax = plot_trajectory_3d(ccar.X, ccar.Y, ccar.Z)
plot_trajectory_3d!(fig, ax, hcar.X, hcar.Y, hcar.Z; linewidth = 3, color=:orange)
fig

hham = compute_invariant(hsol.t, hsol.q, hsol.p, hode.parameters, PauliParticle3d.TokamakIterCylindrical.hamiltonian)
plot_invariant_error(hsol.t, hham; label="H", padding_right=30)

Integrate charged particle dynamics with large time step and convert result to cartesian coordinates:

hode_Δt5 = PauliParticle3d.TokamakIterCylindrical.hodeproblem(X₀, u₀; parameters = parameters, timestep=Δt5, timespan=timespan_long)
hsol_Δt5 = integrate(hode_Δt5, PartitionedGauss(1); initialguess=NoInitialGuess(), options...);
hcar_Δt5 = cartesian_solution(hsol_Δt5, PauliParticle3d.TokamakIterCylindrical);

(t = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0  …  999550.0, 999600.0, 999650.0, 999700.0, 999750.0, 999800.0, 999850.0, 999900.0, 999950.0, 1.0e6], R = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.878511709595285, 6.5772763900007565, 6.240790775295039, 5.968853216207319, 5.801079716553296, 5.726825245819593, 5.741969618975433, 5.845561405663806, 6.05487447244523  …  5.874967190209391, 6.097604184047945, 6.411304887356218, 6.746906299778312, 6.970724097100824, 6.966629130458922, 6.736164717601255, 6.399773075170189, 6.088052134304921, 5.869392497012748], X = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 5.296512003921667, 1.4713924073530837, -2.1207203466298514, -4.261929447704657, -5.14730765004022, -5.391973395995331, -5.364682175129848, -5.033649806537364, -4.056524512533936  …  -5.872718470496306, -5.728947041511643, -4.149260557926953, -0.5671665291807318, 3.9019694117589254, 6.695229191724383, 6.212571263077043, 3.3829468531645297, 0.26676205383568763, -1.8913360085703457], Y = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 4.388722357515127, 6.410582586165245, 5.8694134896425405, 4.1788953204664425, 2.6753971356872417, 1.9295466547962035, 2.0470467178550087, 2.9720292348044905, 4.495121105852458  …  0.1625338497924564, -2.088047552349803, -4.887582958996137, -6.723025118662383, -5.776298914322198, -1.9255719959639253, 2.603819003664678, 5.432565324258651, 6.082204928859934, 5.55631320091612], Z = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.4467052127544391, 0.7495337923642542, 0.8720656301329122, 0.8766972558934716, 0.8475618367712934, 0.8421279965650896, 0.8813005546869146, 0.9598663919873074, 1.0443593713576782  …  -0.861846472211576, -0.8840900054412248, -0.828239888380784, -0.6167362085589778, -0.2256712352756829, 0.24123285578547704, 0.6271892537733403, 0.8323255105541538, 0.8841427978084581, 0.8607389404201203])

Pauli Particle with Variational Integrator

Integrate charged particle dynamics and convert result to cartesian coordinates:

pode = PauliParticle3d.TokamakIterCylindrical.iodeproblem(X₀, u₀; parameters = parameters, timestep=Δt3, timespan=timespan)
psol = integrate(pode, VPRKGauss(1); options...)
pcar = cartesian_solution(psol, PauliParticle3d.TokamakIterCylindrical);

(t = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …  9991.0, 9992.0, 9993.0, 9994.0, 9995.0, 9996.0, 9997.0, 9998.0, 9999.0, 10000.0], R = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 7.000218761141979, 6.999999911762745, 6.99968271654664, 6.999671599626936, 6.999066237995937, 6.998711774229856, 6.998335908170622, 6.997471736575916, 6.997039876571222  …  5.838475038376013, 5.839237543731181, 5.840100721066344, 5.840644002206931, 5.841713597230302, 5.842413315931241, 5.843291051652711, 5.844390660044442, 5.845095810206714, 5.8463212525143495], X = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.99972110610099, 6.998010317076701, 6.995205596450788, 6.991713781674032, 6.986635998983824, 6.980814527465675, 6.973984683829918, 6.965677732307812, 6.956813247227001  …  -5.641146036643115, -5.64378038957999, -5.646598521258924, -5.649194556112123, -5.652367495636728, -5.655274368670241, -5.6584219903997965, -5.661859680360303, -5.66499716509662, -5.668697673733367], Y = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.08346939947349326, 0.16688429155105253, 0.25031339497780813, 0.3336778353042158, 0.41694798423532736, 0.5001949938452154, 0.5833036201636775, 0.6662915531581167, 0.7492105695562088  …  1.5050787909626089, 1.4981448549162697, 1.4908727551056025, 1.4831463945642562, 1.4752489438514367, 1.4668555379559634, 1.4581875369391468, 1.44922294596968, 1.4397750344909195, 1.4301534433097969], Z = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.008225473113949465, 0.016820456442639915, 0.024899349641922502, 0.03331836604043395, 0.0417715511032882, 0.04982403243741868, 0.05838201203211183, 0.06661932922033084, 0.07477324893331222  …  -0.8936857913199071, -0.8935995955901957, -0.8930377620831371, -0.8928316577731438, -0.8925975431841093, -0.8921006127625281, -0.8920896088981191, -0.8916341043420021, -0.8914321801092479, -0.8913281418473725])

Plot solution and energy error:

fig, ax = plot_trajectory_poloidal(ccar.R, ccar.Z, ChargedParticle3d.TokamakIterCylindrical; label="Charged Particle Δt=0.01", linewidth = 1)
plot_trajectory_scatter!(fig, ax, pcar.R, pcar.Z; label="Pauli Particle Δt=1.0")
axislegend(ax)
fig

fig, ax = plot_trajectory_3d(ccar.X, ccar.Y, ccar.Z)
plot_trajectory_3d!(fig, ax, pcar.X, pcar.Y, pcar.Z; linewidth = 3, color=:orange)
fig

pham = compute_invariant(psol.t, psol.q, psol.p, pode.parameters, PauliParticle3d.TokamakIterCylindrical.hamiltonian)
plot_invariant_error(psol.t, pham; label="H", padding_right=30)

Integrate charged particle dynamics with large time step and convert result to cartesian coordinates:

pode_Δt5 = PauliParticle3d.TokamakIterCylindrical.iodeproblem(X₀, u₀; parameters = parameters, timestep=Δt5, timespan=timespan_long)
psol_Δt5 = integrate(pode_Δt5, VPRKGauss(1); initialguess=NoInitialGuess(), options...);
pcar_Δt5 = cartesian_solution(psol_Δt5, PauliParticle3d.TokamakIterCylindrical);

(t = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0  …  999550.0, 999600.0, 999650.0, 999700.0, 999750.0, 999800.0, 999850.0, 999900.0, 999950.0, 1.0e6], R = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.878511709588753, 6.57727639001839, 6.240790775294953, 5.96885321622757, 5.801079716555985, 5.726825245841677, 5.741969618981687, 5.845561405691579, 6.054874472458627  …  5.874967539588218, 6.097604744500593, 6.411305589472413, 6.746906923051515, 6.970724347040928, 6.966628863929134, 6.736164086993547, 6.399772375902961, 6.088051581571603, 5.869392156069396], X = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 5.296512003968445, 1.4713924074516174, -2.1207203465158155, -4.261929447625891, -5.147307649972677, -5.391973395959284, -5.364682175067047, -5.033649806452352, -4.056524512362969  …  -5.8727188768913186, -5.728946180419974, -4.149256527984896, -0.5671593108354429, 3.9019762427837223, 6.695231163024868, 6.212567878251014, 3.382941544961088, 0.2667580592986503, -1.8913377783284593], Y = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 4.388722357448435, 6.4105825861607215, 5.869413489683652, 4.178895320575699, 2.675397135823021, 1.9295466549624776, 2.0470467180371332, 2.9720292350030997, 4.495121106024787  …  0.1625317944883697, -2.088051551568481, -4.887587301171707, -6.723026353098218, -5.776294601488103, -1.9255641774272823, 2.6038254482636316, 5.43256780599625, 6.082204550792066, 5.556312238346268], Z = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.44670521274889907, 0.7495337923614751, 0.872065630133583, 0.8766972558969433, 0.8475618367774139, 0.8421279965729346, 0.8813005546959753, 0.9598663919947872, 1.0443593713625  …  -0.8618465378450402, -0.884090000312148, -0.8282396389294742, -0.616735577865596, -0.22567030553571032, 0.24123378197429995, 0.6271898682603432, 0.8323257511442208, 0.8841427954280274, 0.8607388779751333])

Guiding Center 4D with Variational Lobatto-IIIA-IIIB

The projection is not optional here. The degenerate variational discretisation of the 4D guiding centre carries a parasitic mode that plain VPRK does not control, so VPRKLobattoIIIAIIIB(2) on its own throws on a NaN direction vector at this step, and SymmetricProjection — which constrains exactly the drift off p = ϑ(q) that the mode feeds on — holds the orbit at 3.6E-3 relative energy. SymmetricProjection(VPRKGauss(2)) is better still at 2.5E-7, and is what the test suite uses for every 4D variational orbit; the Lobatto pair is kept here to show what the method itself does. See Findings, "The variational formulation is unstable here, and finer steps make it worse".

Integrate charged particle dynamics and convert result to cartesian coordinates:

vode = GuidingCenter4d.TokamakIterCylindrical.iodeproblem(q₀; parameters = parameters, timestep=10., timespan=(0,1800))
vsol = integrate(vode, SymmetricProjection(VPRKLobattoIIIAIIIB(2)); options...)
vcar = cartesian_solution(vsol, GuidingCenter4d.TokamakIterCylindrical);

(t = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0  …  1710.0, 1720.0, 1730.0, 1740.0, 1750.0, 1760.0, 1770.0, 1780.0, 1790.0, 1800.0], R = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.996119670605115, 6.984555426792706, 6.965534721724515, 6.9394275285841776, 6.906733186157716, 6.86806298426294, 6.824119409553264, 6.775673125604792, 6.72353884303617  …  6.060835360096336, 6.109063539784369, 6.162413710678475, 6.2208169670936115, 6.284107632206666, 6.35200548826252, 6.424099660958921, 6.499835491353432, 6.578505806001775, 6.65924796142018], X = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.94752275204412, 6.791482377597106, 6.535991114852292, 6.18770022106727, 5.75550915944809, 5.250183315014257, 4.683902963088454, 4.069769817678476, 3.421299196345918  …  -3.0424738806770106, -2.7085611021531943, -2.325974727104847, -1.8926786039140053, -1.4077110122479213, -0.8715038056758694, -0.2861892290326275, 0.3441303653512677, 1.0131772416164133, 1.7123849030736433], Y = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.8231755919954897, 1.6308836944207301, 2.408213882968394, 3.141340541659877, 3.8179939261713676, 4.427850981520116, 4.962827696808825, 5.4172612948016745, 5.787977745542819  …  5.241858253290741, 5.475797119052659, 5.706591303959622, 5.925903428201129, 6.1244067826324216, 6.291935698941118, 6.417721728083747, 6.490719197924133, 6.500016193569799, 6.435318279262051], Z = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.082476363371682, 0.16401164537740212, 0.24369223439404633, 0.32065821291307284, 0.39412716761415123, 0.46341456379520085, 0.5279498712064379, 0.587287882784407, 0.6411149500081341  …  1.0277644426964108, 1.0396308989124428, 1.0494607223917052, 1.0565266478530122, 1.0600470404874145, 1.0592005156739421, 1.0531453587231725, 1.0410437298118127, 1.0220902607229079, 0.9955442263713089])

Plot solution and energy error:

fig, ax = plot_trajectory_poloidal(ccar.R, ccar.Z, ChargedParticle3d.TokamakIterCylindrical; label="Charged Particle Δt=0.01", linewidth = 1)
plot_trajectory_scatter!(fig, ax, vcar.R, vcar.Z; label="Guiding Center Δt=10.0", color=:orange, markersize=10)
axislegend(ax)
fig

fig, ax = plot_trajectory_3d(vcar.X, vcar.Y, vcar.Z)
fig

vham = compute_invariant(vsol.t, vsol.q, vsol.p, vode.parameters, GuidingCenter4d.TokamakIterCylindrical.hamiltonian)
plot_invariant_error(vsol.t, vham; label="H", padding_right=30)

Guiding Center 4D with Projected Variational Midpoint

Integrate charged particle dynamics and convert result to cartesian coordinates:

gode = GuidingCenter4d.TokamakIterCylindrical.iodeproblem(q₀; parameters = parameters, timestep=Δt3, timespan=timespan)
gsol = integrate(gode, SymmetricProjection(VPRKGauss(1)); options...)
gcar = cartesian_solution(gsol, GuidingCenter4d.TokamakIterCylindrical);

(t = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …  9991.0, 9992.0, 9993.0, 9994.0, 9995.0, 9996.0, 9997.0, 9998.0, 9999.0, 10000.0], R = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.999960425214316, 6.9998417088443246, 6.999643874848232, 6.999366963146465, 6.9990110296068035, 6.998576146023575, 6.998062400090929, 6.997469895370198, 6.996798751251373  …  5.838639969481633, 5.839360104802827, 5.840116600662457, 5.840909421317787, 5.841738532798407, 5.84260390287886, 5.843505501050233, 5.844443298490703, 5.84541726803506, 5.846427384143189], X = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.9994628713178795, 6.997851631165597, 6.995166717158571, 6.991408858433146, 6.98657907531271, 6.980678678840557, 6.9737092701797865, 6.96567273988062, 6.956571267015607  …  -5.641979458581033, -5.6446031954341045, -5.6473445826393185, -5.650202159880404, -5.6531744065566265, -5.65625974175553, -5.659456524220685, -5.662763052314729, -5.666177563977976, -5.669698236682917], Y = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.08345937699895489, 0.16690265746723962, 0.2503137508355221, 0.3336765784538966, 0.4169750795434686, 0.5001932171381922, 0.5833149840137237, 0.6663244086000604, 0.7492055608747467  …  1.5025925868904615, 1.4955203775468855, 1.4881065399586, 1.4803507761749228, 1.4722527685947167, 1.463812180483191, 1.4550286565082662, 1.4459018232964718, 1.4364312900083607, 1.4266166489334084], Z = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.008335240610268174, 0.01666950807333152, 0.025001829533202822, 0.033331232716193315, 0.04165674622172026, 0.049977399812704285, 0.05829222470542082, 0.0666002538586694, 0.07490052226212561  …  -0.8936953924049514, -0.8933577826690524, -0.8930342488065028, -0.8927246856278694, -0.8924289836192608, -0.8921470289495685, -0.8918787034772548, -0.8916238847566944, -0.8913824460440871, -0.8911542563029528])

Plot solution and energy error:

fig, ax = plot_trajectory_poloidal(ccar.R, ccar.Z, ChargedParticle3d.TokamakIterCylindrical; label="Charged Particle Δt=0.01", linewidth = 1)
plot_trajectory_scatter!(fig, ax, gcar.R, gcar.Z; label="Guiding Center 4D Δt=1.0")
axislegend(ax)
fig

fig, ax = plot_trajectory_3d(ccar.X, ccar.Y, ccar.Z)
plot_trajectory_3d!(fig, ax, gcar.X, gcar.Y, gcar.Z; linewidth = 3, color=:orange)
fig

gham = compute_invariant(gsol.t, gsol.q, gsol.p, gode.parameters, GuidingCenter4d.TokamakIterCylindrical.hamiltonian)
plot_invariant_error(gsol.t, gham; label="H", padding_right=30)

Integrate charged particle dynamics with large time step and convert result to cartesian coordinates:

gode_Δt5 = GuidingCenter4d.TokamakIterCylindrical.iodeproblem(q₀; parameters = parameters, timestep=Δt5, timespan=timespan_long)
gsol_Δt5 = integrate(gode_Δt5, SymmetricProjection(VPRKGauss(1)); options...);
gcar_Δt5 = cartesian_solution(gsol_Δt5, GuidingCenter4d.TokamakIterCylindrical);

(t = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0  …  999550.0, 999600.0, 999650.0, 999700.0, 999750.0, 999800.0, 999850.0, 999900.0, 999950.0, 1.0e6], R = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.877900620346197, 6.578801239133602, 6.240971168280195, 5.969739720406599, 5.800414412755771, 5.727223536751797, 5.741240388861567, 5.846245833439823, 6.054861701253614  …  5.925531936598335, 6.176206373343365, 6.506398763203715, 6.825327180209302, 6.994318788992685, 6.922371506108035, 6.649647403981141, 6.308965905393606, 6.018603249472279, 5.827618644524463], X = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 5.295727034123648, 1.4702354955256822, -2.122103915125231, -4.263611269715783, -5.147271465881278, -5.392690461331111, -5.36421887467768, -5.034501600104523, -4.056570686598178  …  -5.885336026815394, -5.314199766762736, -3.04203520850868, 1.0436637338928525, 5.211512251938284, 6.918144270111571, 5.3733414162776825, 2.150030743608982, -0.7795887181509834, -2.5841403945819934], Y = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 4.388711897973089, 6.412412442421508, 5.869105221131653, 4.178446035328543, 2.6740238997116816, 1.928724715503021, 2.046215303248181, 2.9719327017207977, 4.495062233822819  …  -0.6890200164839984, -3.147190176183524, -5.751456046603537, -6.745061684481522, -4.664829618534127, -0.24188246422353968, 3.91727106815294, 5.931308337707218, 5.96789967284052, 5.22334734505554], Z = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.4446173760028961, 0.7476273632679264, 0.8715089492750715, 0.8765468811889714, 0.8476286576957094, 0.8419720562554536, 0.8814534516205534, 0.9599250761934885, 1.0445995205847463  …  -0.8707915291612836, -0.8794618177398351, -0.7875897541124112, -0.5232749759278784, -0.09922792793883492, 0.3592561781940193, 0.6993590589059552, 0.8582253447222418, 0.8809225655956191, 0.8531217842512155])

Guiding Centre 3D (Modified Poisson Structure)

Integrate charged particle dynamics and convert result to cartesian coordinates:

g3ode = GuidingCenter3d.TokamakIterCylindrical.hodeproblem(q3₀; parameters = parameters, timestep=Δt3, timespan=(0., 19_189.))
g3sol = integrate(g3ode, PartitionedGauss(1); initialguess=NoInitialGuess(), options...)
g3car = cartesian_solution(g3sol, GuidingCenter3d.TokamakIterCylindrical);

(t = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …  19180.0, 19181.0, 19182.0, 19183.0, 19184.0, 19185.0, 19186.0, 19187.0, 19188.0, 19189.0], R = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.9999604252074725, 6.999841708816957, 6.999643874786685, 6.999366963037122, 6.9990110294361045, 6.998576145778032, 6.998062399757142, 6.99746989493487, 6.996798750701324  …  6.511671113760478, 6.517604047135288, 6.523536739381884, 6.5294684642718215, 6.535398489091919, 6.541326074759888, 6.547250475943959, 6.553170941186481, 6.559086713031486, 6.564997028156206], X = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.999462871284684, 6.997851631032899, 6.995166716860321, 6.991408857903721, 6.986579074487087, 6.980678677654478, 6.973709268569928, 6.9656727377847565, 6.956571264372774  …  -6.484288303073909, -6.483791635114939, -6.482590680144662, -6.480677851945411, -6.478045684619624, -6.474686836540002, -6.47059409429192, -6.465760376605735, -6.460178738276533, -6.453842374068916], Y = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.08345937920898217, 0.1669026618831998, 0.25031375744923506, 0.3336765872531116, 0.4169750905118865, 0.5001932302554973, 0.5833149992556224, 0.6663244259383252, 0.7492055802772712  …  -0.5965449659510256, -0.6630298241013226, -0.7297600042252089, -0.7967264306049066, -0.8639198569346574, -0.9313308676669269, -0.9989498794548834, -1.0667671426914958, -1.1347727431455172, -1.2029566036945465], Z = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.008335241782371447, 0.016669510417032594, 0.02500183304749136, 0.03333123739955533, 0.04165675207213978, 0.049977406827665743, 0.058292232881912015, 0.0666002631931849, 0.07490053275067088  …  -0.7906542533219142, -0.787590695530333, -0.7844726174228518, -0.7812997321822179, -0.7780717621259111, -0.7747884388720143, -0.7714495035036046, -0.7680547067315678, -0.7646038090557392, -0.7610965809242735])

Plot solution and energy error:

fig, ax = plot_trajectory_poloidal(ccar.R, ccar.Z, ChargedParticle3d.TokamakIterCylindrical; label="Charged Particle Δt=0.01", linewidth = 1)
plot_trajectory_scatter!(fig, ax, g3car.R, g3car.Z; label="Guiding Center 3D Δt=1.0", color=:orange)
axislegend(ax)
fig

fig, ax = plot_trajectory_3d(g3car.X, g3car.Y, g3car.Z)
fig

g3ham = compute_invariant(g3sol.t, g3sol.q, g3sol.p, g3ode.parameters, GuidingCenter3d.TokamakIterCylindrical.hamiltonian)
plot_invariant_error(g3sol.t, g3ham; label="H", padding_right=40)

fig = plot_invariant_error(g3sol.t, g3ham; label="H", padding_right=40)
ylims!(-0.1, 0.2)
fig

Compute and plot constraints:

g3g1 = compute_invariant(g3sol.t, g3sol.q, g3sol.p, g3ode.parameters, GuidingCenter3d.TokamakIterCylindrical.g₁)
println()
println("g₁(0) = ", g3g1[begin])
println("g₁(T) = ", g3g1[end])
println()
plot_invariant(g3sol.t, g3g1; label="g₁", padding_right=40)

g3g2 = compute_invariant(g3sol.t, g3sol.q, g3sol.p, g3ode.parameters, GuidingCenter3d.TokamakIterCylindrical.g₂)
println()
println("g₂(0) = ", g3g2[begin])
println("g₂(T) = ", g3g2[end])
println()
plot_invariant(g3sol.t, g3g2; label="g₂", padding_right=40)

Integrate charged particle dynamics with large time step and convert result to cartesian coordinates:

g3ode_Δt4 = GuidingCenter3d.TokamakIterCylindrical.hodeproblem(q3₀; parameters = parameters, timestep=Δt4, timespan=(0, 100*Δt4))
g3sol_Δt4 = integrate(g3ode_Δt4, PartitionedGauss(1); initialguess=MidpointExtrapolation(5, 1.0), options...)
g3car_Δt4 = cartesian_solution(g3sol_Δt4, GuidingCenter3d.TokamakIterCylindrical);

(t = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0  …  910.0, 920.0, 930.0, 940.0, 950.0, 960.0, 970.0, 980.0, 990.0, 1000.0], R = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.996013141220926, 6.984133314228864, 6.9645997511586835, 6.937801291376988, 6.904262149855647, 6.864623127106668, 6.819619293659746, 6.770055342582221, 6.716779879620805  …  5.88810747781637, 5.867961670392674, 5.852087917528888, 5.84033968287838, 5.832584332477335, 5.82870685559188, 5.828612129178678, 5.832225853175634, 5.839494263464332, 5.850382704553622], X = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.945837506345509, 6.784827421067112, 6.521334070687655, 6.162409453630437, 5.7174752154656465, 5.197905822066907, 4.6165386839059614, 3.987140798497575, 3.323863050163772  …  -4.960608584067725, -5.038783264096781, -5.0972642468098615, -5.138528673286372, -5.164443565681982, -5.17624752956102, -5.174545407339151, -5.159313142341956, -5.129911284461739, -5.085106734754262], Y = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.8363858006804605, 1.6565732755472231, 2.4449645462312883, 3.1871298192007202, 3.8703117438347756, 4.483840512566285, 5.01943999769457, 5.471412760392164, 5.836699955885591  …  -3.1720927076601417, -3.00726410257825, -2.8748617693171763, -2.775804511347607, -2.710638790459929, -2.6796054039983614, -2.6826850690668302, -2.71962245572944, -2.78992893575893, -2.892864892437426], Z = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.0835548378763297, 0.16613090685472598, 0.24677874386860613, 0.3246061196313822, 0.39880330143390613, 0.4686645373330689, 0.5336048911112099, 0.5931718496105876, 0.6470514423535675  …  -0.9580808915097043, -0.9449465234738347, -0.9328977628334993, -0.9221278335486516, -0.9127741127954397, -0.9049215380936153, -0.8986055822303245, -0.8938146813210838, -0.8904920674106971, -0.8885370160130576])

g3ode_Δt5 = GuidingCenter3d.TokamakIterCylindrical.hodeproblem(q3₀; parameters = parameters, timestep=Δt5, timespan=(0, 100*Δt4))
g3sol_Δt5 = integrate(g3ode_Δt5, PartitionedGauss(1); initialguess=MidpointExtrapolation(5, 1.0), options...)
g3car_Δt5 = cartesian_solution(g3sol_Δt5, GuidingCenter3d.TokamakIterCylindrical);

(t = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0  …  550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0, 900.0, 950.0, 1000.0], R = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 6.877833682436445, 6.578774968520926, 6.24131825165556, 5.970515719874683, 5.801535349148304, 5.728715156534326, 5.743255036115957, 5.84899382883064, 6.058462772205881  …  6.759518914841709, 7.122229405249388, 7.339503726102351, 7.325781387228392, 7.086670181460661, 6.715304416271004, 6.33394875756178, 6.0292472652663935, 5.831841065809858, 5.736839934141012], X = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[7.0, 5.295473040163581, 1.4719688115486231, -2.115969623551643, -4.25589107793536, -5.140674816393608, -5.3869632054706695, -5.356707966819458, -5.020621319302854, -4.0310767551689395  …  1.2737337274093805, 4.805421169202781, 7.1003302545158675, 6.944662153952399, 4.4236224042020185, 0.8517981624886704, -2.265210546222449, -4.204637057827037, -5.092722612131737, -5.378934284095873], Y = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 4.388913469659427, 6.411987827831993, 5.8716885195534285, 4.18741556261388, 2.6891029804590767, 1.9490522742120788, 2.071390394900066, 3.0006816522018096, 4.522763685639703  …  6.6384258941229914, 5.2568126357686396, 1.8583931829074636, -2.332110739421083, -5.536466299531285, -6.661062474828605, -5.91504252264521, -4.321209298063435, -2.8415746713307333, -1.994592288497345], Z = GeometricSolutions.ScalarDataSeries{Float64} with data type Float64
[0.0, 0.4448263312000822, 0.7479237250679251, 0.8718812280821506, 0.8771942603710232, 0.8486521984571659, 0.8433499158300853, 0.8830672326622178, 0.961501480944872, 1.0456396550624292  …  0.9712113885741062, 0.6827521063109769, 0.22920220645594125, -0.28569818415808934, -0.7246723389860626, -0.9917973524578181, -1.0767314055463182, -1.0378691782325538, -0.951716347524209, -0.8764029293134865])

Plot solution, energy error and constraints:

fig, ax = plot_trajectory_poloidal(ccar.R, ccar.Z, ChargedParticle3d.TokamakIterCylindrical; label="Charged Particle Δt=0.01", linewidth = 1)
plot_trajectory_scatter!(fig, ax, g3car.R, g3car.Z; label="Guiding Center 3D Δt=1.0", markersize=12, color=Makie.wong_colors()[6])
plot_trajectory_scatter!(fig, ax, g3car_Δt4.R, g3car_Δt4.Z; label="Guiding Center 3D Δt=10.0", markersize=12, color=Makie.wong_colors()[3])
plot_trajectory_scatter!(fig, ax, g3car_Δt5.R, g3car_Δt5.Z; label="Guiding Center 3D Δt=50.0", markersize=12, color=Makie.wong_colors()[2])
axislegend(ax)
fig

g3ham_Δt4 = compute_invariant(g3sol_Δt4.t, g3sol_Δt4.q, g3sol_Δt4.p, g3ode_Δt4.parameters, GuidingCenter3d.TokamakIterCylindrical.hamiltonian)
plot_invariant_error(g3sol_Δt4.t, g3ham_Δt4; label="H", padding_right=40)

g3g1_Δt4 = compute_invariant(g3sol_Δt4.t, g3sol_Δt4.q, g3sol_Δt4.p, g3ode_Δt4.parameters, GuidingCenter3d.TokamakIterCylindrical.g₁)
plot_invariant(g3sol_Δt4.t, g3g1_Δt4; label="g₁", padding_right=40)

g3g2_Δt4 = compute_invariant(g3sol_Δt4.t, g3sol_Δt4.q, g3sol_Δt4.p, g3ode_Δt4.parameters, GuidingCenter3d.TokamakIterCylindrical.g₂)
plot_invariant(g3sol_Δt4.t, g3g2_Δt4; label="g₂", padding_right=40)

g3ham_Δt5 = compute_invariant(g3sol_Δt5.t, g3sol_Δt5.q, g3sol_Δt5.p, g3ode_Δt5.parameters, GuidingCenter3d.TokamakIterCylindrical.hamiltonian)
plot_invariant_error(g3sol_Δt5.t, g3ham_Δt5; label="H", padding_right=40)

g3g1_Δt5 = compute_invariant(g3sol_Δt5.t, g3sol_Δt5.q, g3sol_Δt5.p, g3ode_Δt5.parameters, GuidingCenter3d.TokamakIterCylindrical.g₁)
plot_invariant(g3sol_Δt5.t, g3g1_Δt5; label="g₁", padding_right=40)

g3g2_Δt5 = compute_invariant(g3sol_Δt5.t, g3sol_Δt5.q, g3sol_Δt5.p, g3ode_Δt5.parameters, GuidingCenter3d.TokamakIterCylindrical.g₂)
plot_invariant(g3sol_Δt5.t, g3g2_Δt5; label="g₂", padding_right=40)