bondscell_results$081ce252-d083-4f03-a5bc-a4b171d222bcqueued¤logsrunning¦outputbody]

Procedures to obtain ion surface flux

  1. Obtain magnetic field, density, velocity and current from the unstructured grid.

Requires a conversion from BATSRUS unstructured *.out format to VTK *.vtu format.

  1. Interpolate MHD quantities onto a Cartesian grid and calculate the derived electric field including the convective term and Hall term.

$$\mathbf{E} = -\sum_\alpha\mathbf{U}_\alpha\times\mathbf{B} + \frac{\mathbf{J}\times\mathbf{B}}{\sum_\alpha q_\alpha n_\alpha} - \frac{\nabla P_e}{\sum_\alpha q_\alpha n_\alpha} + \eta \mathbf{J}$$

This step is performed in ParaView with Python.

It would be better to first convert the basic quantities from spherical coordinates to Cartesian coordinates, then calculate the derived electric field. If we first calculate the electric field and then interpolate to Cartesian coordinates, then around the inner boundary region there might be some artifacts due to interpolations.

  1. Read the field info into TestParticle.jl.

  2. Tracing with prescribed seeds.

Potential Issues

mimetext/htmlrootassigneelast_run_timestampA&~hpersist_js_state·has_pluto_hook_features§cell_id$081ce252-d083-4f03-a5bc-a4b171d222bcdepends_on_disabled_cells§runtimeGpublished_object_keysdepends_on_skipped_cells§errored$e73e1a00-3706-42e1-acd1-ebc8b2b95c62queued¤logsrunning¦outputbodyrows0.0text/plain1.0text/plain0.0text/plain0.0text/plain0.0text/plain1.0text/plain0.1text/plain0.00124295text/plain1.0text/plain0.00124295text/plain0.000124295text/plain0.00118993text/plain0.999999text/plain0.1text/plain0.00902976text/plain1.00004text/plain0.00902965text/plain0.000902976text/plain0.00864452text/plain0.999963text/plain0.1text/plain0.0527251text/plain1.00133text/plain0.0527027text/plain0.00527251text/plain0.0504548text/plain0.998726text/plain0.1text/plain0.301738text/plain1.04328text/plain0.297559text/plain0.0301738text/plain0.284868text/plain0.958567text/plain0.1text/plain1.05209text/plain1.48653text/plain0.883008text/plain0.105209text/plain0.845346text/plain0.534219text/plain0.1text/plain2.14553text/plain2.52988text/plain0.924956text/plain0.214553text/plain0.885504text/plain-0.464631text/plain0.1text/plain3.30351text/plain3.08887text/plain-0.0219528text/plain0.330351text/plain-0.0210164text/plain-0.999779text/plain0.1text/plain 4.58819text/plain2.37303text/plain-0.99156text/plain0.458819text/plain-0.949268text/plain-0.314468text/plain0.1text/plain 5.9665text/plain1.16581text/plain-0.564705text/plain0.59665text/plain-0.540619text/plain0.841267text/plain0.1text/plainmoreobjectid872690c148ccb9abschemanamestimestampvalue1value2value3value4value5value6typesFloat64Float64Float64Float64Float64Float64Float64mime"application/vnd.pluto.table+objectrootassigneelast_run_timestampA&persist_js_state·has_pluto_hook_features§cell_id$e73e1a00-3706-42e1-acd1-ebc8b2b95c62depends_on_disabled_cells§runtimexpublished_object_keysdepends_on_skipped_cells§errored$d8d1dc93-44e7-466f-93fe-a9b427368944queued¤logsrunning¦outputbody

WIP: Tracing Particles Near Ganymede

Ganymede is believed to have both an ionosphere and an exosphere besides the magnetosphere.

Our ultimate goal is to build a self-consistent magnetosphere-ionosphere-exosphere coupling model for Ganymede.

mimetext/htmlrootassigneelast_run_timestampA&qpersist_js_state·has_pluto_hook_features§cell_id$d8d1dc93-44e7-466f-93fe-a9b427368944depends_on_disabled_cells§runtime!published_object_keysdepends_on_skipped_cells§errored$9c31a2c5-e0c7-4891-ad45-7a0e10012bd8queued¤logsrunning¦outputbody

Relativistic B-only tracing:

$$\begin{aligned} \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} &= \frac{q}{\gamma m c}\mathbf{u}\times\mathbf{B} \\ \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} &= \mathbf{u} \end{aligned}$$

Initial particle seeds: random location, $u = c$, pitch angle $\theta = \cos^{-1}(\mu_0),\,\text{where }\mu_0 = 0.1-0.7$, random azimuthal angle.

The Lorentz equation in SI units is written as

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t} = \frac{q}{m}\left( \mathbf{v}\times\mathbf{B} + \mathbf{E} \right)$$

It can be normalized to

$$\frac{\mathrm{d}\mathbf{v}^\prime}{\mathrm{d}t^\prime} = \mathbf{v}^\prime\times\mathbf{B}^\prime + \mathbf{E}^\prime$$

with the following transformation

$$\begin{aligned} \mathbf{v} &= \mathbf{v}^\prime V_0 \\ \mathbf{x} &= \mathbf{x}^\prime l_0 \\ t &= t^\prime t_0 \\ \mathbf{B} &= \mathbf{B}^\prime B_0 \\ \mathbf{E} &= \mathbf{E}^\prime E_0 = \mathbf{E}^\prime V_0 B_0 \end{aligned}$$

where $V_0, l_0, t_0, B_0$ are prescribed parameters.

After the normalization, we can think of everything in the dimensionless natural units. If the magnetic field is homogeneous and the initial perpendicular velocity is 1, then the gyroradius is 1. MHD solutions are also dimensionless by nature. For the turbulence EM fields, we know their values and the spatial range. To make them more easily understandable, we can convert from the dimensionless MHD to natural units via

$$\begin{aligned} V_0 &= c \\ l_0 &= 4 / nx \\ t_0 &= l_0 / V_0 \\ B_0 &= \mathrm{norm}(B) \end{aligned}$$

and

x = range(xmin, xmax, length=nx) / l₀ 
y = range(ymin, ymax, length=ny) / l₀
z = range(zmin, zmax, length=nz) / l₀
B ./= B₀

When the X extent is 1, then after the normalization, the x grid resolution becomes 1/4, which is a quarter of the gyroradius.

mimetext/htmlrootassigneelast_run_timestampA&~h; persist_js_state·has_pluto_hook_features§cell_id$9c31a2c5-e0c7-4891-ad45-7a0e10012bd8depends_on_disabled_cells§runtime ҵpublished_object_keysdepends_on_skipped_cells§errored$9643aa92-18fe-4850-a246-8ad134c085e8queued¤logsrunning¦outputbodyY

The surface flux is very sensitive to the seeding procedures. Placing a plane source in the upstream as in [Plainaki+ 2015, 2020] is deficient in capturing all the possible sources, which ends up in significantly underestimating ion sputtering.

Stop Conditions

mimetext/htmlrootassigneelast_run_timestampA&~h2ʰpersist_js_state·has_pluto_hook_features§cell_id$9643aa92-18fe-4850-a246-8ad134c085e8depends_on_disabled_cells§runtimepublished_object_keysdepends_on_skipped_cells§errored$45a18eea-80c0-4aa4-86d4-c25f90df9303queued¤logsrunning¦outputbody;

Julia Community Ecosystem

Advantages of Julia's Open Source Community

mimetext/htmlrootassigneelast_run_timestampA&~h~:persist_js_state·has_pluto_hook_features§cell_id$45a18eea-80c0-4aa4-86d4-c25f90df9303depends_on_disabled_cells§runtimeYlpublished_object_keysdepends_on_skipped_cells§errored$85c0b164-b6e0-43f5-a458-939b4fd7dddequeued¤logsrunning¦outputbody

$$\mu=0.1,\,l_0=16/nx$$

mimetext/htmlrootassigneelast_run_timestampA&Xְpersist_js_state·has_pluto_hook_features§cell_id$85c0b164-b6e0-43f5-a458-939b4fd7dddedepends_on_disabled_cells§runtimepublished_object_keysdepends_on_skipped_cells§errored$b7244973-218a-4bb6-8a8f-95f84a8fb225queued¤logsrunning¦outputbody

Magnetospheric Ion Motion

mimetext/htmlrootassigneelast_run_timestampA&&ٰpersist_js_state·has_pluto_hook_features§cell_id$b7244973-218a-4bb6-8a8f-95f84a8fb225depends_on_disabled_cells§runtimeXpublished_object_keysdepends_on_skipped_cells§errored$d058eddf-a8b1-453b-a440-885a3da0a1a0queued¤logsrunning¦outputbodyL mimetext/htmlrootassigneelast_run_timestampA&t"persist_js_state·has_pluto_hook_features§cell_id$d058eddf-a8b1-453b-a440-885a3da0a1a0depends_on_disabled_cells§runtime͈published_object_keysdepends_on_skipped_cells§errored$72412346-af0c-49dc-8c0d-c5b4fe97429bqueued¤logsrunning¦outputbody

Our Model

We use ideal MHD G2 spherical output as a starting point.

Seeding

An important question is how to set seeds.

mimetext/htmlrootassigneelast_run_timestampA&]ٰpersist_js_state·has_pluto_hook_features§cell_id$72412346-af0c-49dc-8c0d-c5b4fe97429bdepends_on_disabled_cells§runtime published_object_keysdepends_on_skipped_cells§errored$bf6d18d0-ee99-4a2e-badb-c1f42071903bqueued¤logslinemsg;Precompiling TestParticle... 437.0 ms ✓ JLLWrappers 523.9 ms ✓ PrecompileTools 660.0 ms ✓ TruncatedStacktraces 700.1 ms ✓ OpenSpecFun_jll 704.5 ms ✓ Rmath_jll 1113.9 ms ✓ Rmath 2161.2 ms ✓ RecipesBase 2855.1 ms ✓ SpecialFunctions 1262.5 ms ✓ DualNumbers 2834.0 ms ✓ RecursiveArrayTools 2277.3 ms ✓ SpecialFunctions → SpecialFunctionsChainRulesCoreExt 1285.0 ms ✓ HypergeometricFunctions 1616.0 ms ✓ StatsFuns 2346.2 ms ✓ StatsFuns → StatsFunsChainRulesCoreExt 10176.9 ms ✓ StaticArrays 834.7 ms ✓ StaticArrays → StaticArraysStatisticsExt 1013.2 ms ✓ StaticArrays → StaticArraysChainRulesCoreExt 1108.6 ms ✓ ConstructionBase → ConstructionBaseStaticArraysExt 862.7 ms ✓ Adapt → AdaptStaticArraysExt 8647.3 ms ✓ SciMLBase 3359.2 ms ✓ Rotations 3369.9 ms ✓ NearestNeighbors 3289.6 ms ✓ Interpolations 1528.8 ms ✓ SciMLBase → SciMLBaseChainRulesCoreExt 7832.1 ms ✓ Distributions 861.9 ms ✓ Distributions → DistributionsChainRulesCoreExt 2864.4 ms ✓ Meshes 7147.0 ms ✓ TestParticle 28 dependencies successfully precompiled in 26 seconds. 103 already precompiled. 1 dependency precompiled but a different version is currently loaded. Restart julia to access the new version. Otherwise, loading dependents of this package may trigger further precompilation to work with the unexpected version. Precompiling OrdinaryDiffEq... 880.8 ms ✓ IntelOpenMP_jll 1461.7 ms ✓ FiniteDiff → FiniteDiffStaticArraysExt 1721.0 ms ✓ CPUSummary 2236.0 ms ✓ StaticArrayInterface 1780.0 ms ✓ ArnoldiMethod 962.1 ms ✓ PolyesterWeave 2270.1 ms ✓ LazyArrays → LazyArraysStaticArraysExt 1965.1 ms ✓ MKL_jll 1247.5 ms ✓ StaticArrayInterface → StaticArrayInterfaceOffsetArraysExt 889.7 ms ✓ CloseOpenIntervals 1497.4 ms ✓ StaticArrayInterface → StaticArrayInterfaceStaticArraysExt 1024.3 ms ✓ LayoutPointers 5986.0 ms ✓ ForwardDiff 1917.0 ms ✓ StrideArraysCore 1373.7 ms ✓ ForwardDiff → ForwardDiffStaticArraysExt 1399.3 ms ✓ PreallocationTools 1210.4 ms ✓ Polyester 1683.3 ms ✓ NLSolversBase 1109.3 ms ✓ FastBroadcast 6642.5 ms ✓ Graphs 1318.0 ms ✓ RecursiveArrayTools → RecursiveArrayToolsFastBroadcastExt 1450.1 ms ✓ VertexSafeGraphs 3005.6 ms ✓ LineSearches 2768.7 ms ✓ SparseDiffTools 12544.6 ms ✓ VectorizationBase 1490.1 ms ✓ SLEEFPirates 7681.6 ms ✓ DiffEqBase 43693.8 ms ✓ LoopVectorization 1916.9 ms ✓ LoopVectorization → SpecialFunctionsExt 2358.0 ms ✓ TriangularSolve 48781.8 ms ✓ SimpleNonlinearSolve 2325.2 ms ✓ SimpleNonlinearSolve → SimpleNonlinearSolveStaticArraysExt 9235.7 ms ✓ RecursiveFactorization 76765.1 ms ✓ ExponentialUtilities 34376.6 ms ✓ LinearSolve 1935.6 ms ✓ LinearSolve → LinearSolveRecursiveArrayToolsExt 90560.9 ms ✓ NonlinearSolve 74759.3 ms ✓ OrdinaryDiffEq 38 dependencies successfully precompiled in 277 seconds. 125 already precompiled. Precompiling DiffEqBaseUnitfulExt... 1167.5 ms ✓ DiffEqBase → DiffEqBaseUnitfulExt 1 dependency successfully precompiled in 1 seconds. 111 already precompiled. Precompiling DiffEqBaseChainRulesCoreExt... 1052.1 ms ✓ DiffEqBase → DiffEqBaseChainRulesCoreExt 1 dependency successfully precompiled in 1 seconds. 114 already precompiled. Precompiling DiffEqBaseDistributionsExt... 877.9 ms ✓ Distributions → DistributionsTestExt 1109.1 ms ✓ DiffEqBase → DiffEqBaseDistributionsExt 2 dependencies successfully precompiled in 1 seconds. 123 already precompiled. Precompiling ForwardDiffExt... 1062.1 ms ✓ LoopVectorization → ForwardDiffExt 1 dependency successfully precompiled in 1 seconds. 77 already precompiled. text/plaincell_id$bf6d18d0-ee99-4a2e-badb-c1f42071903bkwargsidPlutoRunner_d1acb81efileP/home/runner/.julia/packages/Pluto/u5fFh/src/runner/PlutoRunner/src/io/stdout.jlgroupstdoutlevelLogLevel(-555)running¦outputbodymimetext/plainrootassigneelast_run_timestampA&kpersist_js_state·has_pluto_hook_features§cell_id$bf6d18d0-ee99-4a2e-badb-c1f42071903bdepends_on_disabled_cells§runtimeH@Xpublished_object_keysdepends_on_skipped_cells§errored$05aee60b-fdb4-4856-a1cd-a0dbb470bbb3queued¤logsrunning¦outputbody

Preliminary Results

mimetext/htmlrootassigneelast_run_timestampA&꠲persist_js_state·has_pluto_hook_features§cell_id$05aee60b-fdb4-4856-a1cd-a0dbb470bbb3depends_on_disabled_cells§runtime ="published_object_keysdepends_on_skipped_cells§errored$36be1db4-39d5-4067-b42d-53b669bf4ae3queued¤logsrunning¦outputbody_

Motion of Charged Particles

Non-relativistic case

$$m\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t} = q(\mathbf{E} + \mathbf{v}\times\mathbf{B}) \tag{1}$$

Relativistic case

The four-dimensional momentum equation satisfying Lorentz covariance is

$$\frac{\mathrm{d}p^\mu}{\mathrm{d}\tau} = \gamma qF^{\mu \nu}v_\nu \tag{2}$$

where $\gamma = 1/\sqrt{1-v^2/c^2}$ is the Lorentz factor， $\tau$ is the proper time， $p^\mu = \gamma m v^\mu$ is the four-dimensional momentum, and $F^{\mu\nu}$ is the electromagnetic tensor。

There are two forms that can be used in numerical computation after simplification:

$$\frac{\mathrm{d}(\gamma v^i)}{\mathrm{d}t} = \frac{q}{m}(E^i+\epsilon^{ijk}v_j B_k) \tag{3}$$

or

$$\frac{\mathrm{d}v^i}{\mathrm{d}t} = \frac{q}{m \gamma^3}(E^i+\epsilon^{ijk}v_j B_k) \tag{4}$$

Currently we only include Eq.(4) in TestParticle.jl. This needs further testing.

mimetext/htmlrootassigneelast_run_timestampA&~h}8persist_js_state·has_pluto_hook_features§cell_id$36be1db4-39d5-4067-b42d-53b669bf4ae3depends_on_disabled_cells§runtimeõpublished_object_keysdepends_on_skipped_cells§errored$d50e4350-49ac-4203-870b-903a8801156equeued¤logsrunning¦outputbody,mimetext/htmlrootassigneelast_run_timestampA&~h|persist_js_state·has_pluto_hook_features§cell_id$d50e4350-49ac-4203-870b-903a8801156edepends_on_disabled_cells§runtimegpublished_object_keysdepends_on_skipped_cells§errored$45b02635-24ef-4a95-adbd-f6fee4413f21queued¤logsrunning¦outputbody

$$\mu=0.7,\,l_0=4/nx$$

mimetext/htmlrootassigneelast_run_timestampA&|persist_js_state·has_pluto_hook_features§cell_id$45b02635-24ef-4a95-adbd-f6fee4413f21depends_on_disabled_cells§runtimepublished_object_keysdepends_on_skipped_cells§errored$84dcadfc-4b85-48fc-b7a5-5bbc71d66b7fqueued¤logsrunning¦outputbody/

Demos

Electron and Ion Motion

Grad-B and Curvature Drift

Finite-Larmor-Radius Effect

Polarization Drift

Adiabatic motion

mimetext/htmlrootassigneelast_run_timestampA&ܝpersist_js_state·has_pluto_hook_features§cell_id$84dcadfc-4b85-48fc-b7a5-5bbc71d66b7fdepends_on_disabled_cells§runtime 3&published_object_keysdepends_on_skipped_cells§errored$75936fde-d5a7-4920-a4ca-a626a2e90123queued¤logsrunning¦outputbody9

What Are Test Particles?

A test particle model describes non-self-consistent field-particle interactions. A complete description of the particle phase space trajectory not only needs to reflect the effect of the field on the particle, but also needs to reflect the effect of the particle on the field.

Field <–> Particle

However, as an approximate treatment, we can ignore the feedback of a small number of particles on the field and only focus on the effect of the field on the particle.

Field –> Particle

Applications of Test Particle Models

mimetext/htmlrootassigneelast_run_timestampA&Opersist_js_state·has_pluto_hook_features§cell_id$75936fde-d5a7-4920-a4ca-a626a2e90123depends_on_disabled_cells§runtimeZpublished_object_keysdepends_on_skipped_cells§errored$0906fe4a-8f2d-4e6d-8726-97bc493da577queued¤logsrunning¦outputbodyl

Scales in the System

$$\mu=0.1,\,l_0=4/nx$$

mimetext/htmlrootassigneelast_run_timestampA&3persist_js_state·has_pluto_hook_features§cell_id$0906fe4a-8f2d-4e6d-8726-97bc493da577depends_on_disabled_cells§runtimepublished_object_keysdepends_on_skipped_cells§errored$db26a01b-b01c-41dc-b23a-02dbb3867655queued¤logsrunning¦outputbodyK

Demo: Tracing Electrons in Mars' Crustal Field

Goal: Given Mars' crustal field and the locations of the observed electrons, can we determine their source location?

Unknown factors:

mimetext/htmlrootassigneelast_run_timestampA&mpersist_js_state·has_pluto_hook_features§cell_id$db26a01b-b01c-41dc-b23a-02dbb3867655depends_on_disabled_cells§runtime j*published_object_keysdepends_on_skipped_cells§errored$47fd7e19-9859-404e-a0e0-be356eef83b0queued¤logsrunning¦outputbody

$$\mu=0.7,\,l_0=16/nx$$

mimetext/htmlrootassigneelast_run_timestampA&#persist_js_state·has_pluto_hook_features§cell_id$47fd7e19-9859-404e-a0e0-be356eef83b0depends_on_disabled_cells§runtimeepublished_object_keysdepends_on_skipped_cells§errored$b3e0e1b5-4f11-41b4-b2e3-fdb36c6c8b20queued¤logsrunning¦outputbody

Demo: Tracing Particles in MHD Turbulence

$$\begin{aligned} B_\mathrm{in-plane} &= 1 \\ B_z &= 4 \\ U_\mathrm{in-plane} &= 1 \\ U_z &= 0 \\ \end{aligned}$$

mimetext/htmlrootassigneelast_run_timestampA&ϰpersist_js_state·has_pluto_hook_features§cell_id$b3e0e1b5-4f11-41b4-b2e3-fdb36c6c8b20depends_on_disabled_cells§runtime1published_object_keysdepends_on_skipped_cells§errored$96a6b000-b558-11ee-1024-4d2ef81f5ad0queued¤logsrunning¦outputbody٨

TestParticle.jl: A New Tool for An Old Problem

2024/01/31 Hongyang Zhou

mimetext/htmlrootassigneelast_run_timestampA&~h}persist_js_state·has_pluto_hook_features§cell_id$96a6b000-b558-11ee-1024-4d2ef81f5ad0depends_on_disabled_cells§runtimepublished_object_keysdepends_on_skipped_cells§errored$021d3bb6-ceb6-460a-a56e-c8b094f3d932queued¤logsrunning¦outputbody

The Goal of TestParticle.jl

  1. Simplicity: Users only need to provide field information and particle initial conditions

  2. Generality: Compatible with general models, can run independently

  3. Efficiency: The efficiency of C

  4. Visuality: Connected to powerful visualization tools

mimetext/htmlrootassigneelast_run_timestampA&~h~kpersist_js_state·has_pluto_hook_features§cell_id$021d3bb6-ceb6-460a-a56e-c8b094f3d932depends_on_disabled_cells§runtimepublished_object_keysdepends_on_skipped_cells§errored$e8daa564-73a1-4306-87ab-7334cb44f96dqueued¤logsrunning¦outputbody

Design

DifferentialEquations.jl Interface

param = prepare(E, B, species)
prob = ODEProblem(trace!, stateinit, tspan, param)
sol = solve(prob, algorithm)

By taking advantage of Julia's multiple dispatch property, we come up with 4 ODE functions:

Energy Conservation

We need to make sure that the energy is conserved in a B-only tracing for validity. In classical mechanics, Hamiltonian's equations naturally have that the solutions reside on a symplectic manifold in phase space, with the natural splitting being the position and momentum components. For the true Hamiltonian solution, that phase space path is constant energy.

For example, the most commonly used solver in DifferentialEquations.jl, Tsit5() (adaptive 4/5th order Runge-Kutta method), does not conserve energy well. The following shows the test result of a proton in a static magnetic field gyrating for 304 cycles. All adaptive timestepping schemes from DifferentialEquations.jl use a relative tolerence of 1e-3 and an absolute tolerence of 1e-6.

SchemeEnergy difference ratioTime [ms]Memory [MiB]fixed time step
Boris0%0.20.047true
ImplicitMidpoint0%6.12.974true
SSPSDIRK20%11.74.090true
Vern80.17%0.91.707false
RosenbrockW6S4OS0.23%11.53.464true
Vern9-0.59%1.11.662false
Trapezoid1.14%20.25.088false
Vern6-5.36%1.12.036false
Tsit520.59%1.22.309false
BS3-89.73%2.63.503false

Native Boris Solver

The classical Boris method is not a general ODE solver, but designed specifically for the Lorentz equation with the symmetry of the EM fields.

param = prepare(E, B, species=Electron)
prob = TraceProblem(stateinit, tspan, param)
sol = trace_trajectory(prob; dt, savestepinterval)

This method is 2nd order accurate in space. A fixed dt is required to conserve phase space volume, even though it is not a symplectic method (Qin+ 2013). It is the de facto algorithm for particle pusher after being proposed in the 1960s. However, the phase error may still be an issue.

Implementation: Analytic or Numerical Fields?

function getinterp(A, gridx, gridy, gridz, order::Int=1, bc::Int=1)
   @assert size(A,1) == 3 && ndims(A) == 4 "Inconsistent 3D force field and grid!"

   Ax = @view A[1,:,:,:]
   Ay = @view A[2,:,:,:]
   Az = @view A[3,:,:,:]

   itpx, itpy, itpz = _getinterp(Ax, Ay, Az, order, bc)

   interpx = scale(itpx, gridx, gridy, gridz)
   interpy = scale(itpy, gridx, gridy, gridz)
   interpz = scale(itpz, gridx, gridy, gridz)

   # Return field value at a given location.
   function get_field(xu)
      r = @view xu[1:3]

      return SA[interpx(r...), interpy(r...), interpz(r...)]
   end

   return get_field
end
mimetext/htmlrootassigneelast_run_timestampA&~hpersist_js_state·has_pluto_hook_features§cell_id$e8daa564-73a1-4306-87ab-7334cb44f96ddepends_on_disabled_cells§runtime0bhpublished_object_keysdepends_on_skipped_cells§errored±cell_dependencies$081ce252-d083-4f03-a5bc-a4b171d222bcprecedence_heuristic cell_id$081ce252-d083-4f03-a5bc-a4b171d222bcdownstream_cells_mapupstream_cells_map@md_strgetindex$e73e1a00-3706-42e1-acd1-ebc8b2b95c62precedence_heuristic cell_id$e73e1a00-3706-42e1-acd1-ebc8b2b95c62downstream_cells_mapBparamx0tspanstateinitsolEprobv0upstream_cells_mapProtonODEProblemprepareVern9solveSAtrace!$d8d1dc93-44e7-466f-93fe-a9b427368944precedence_heuristic cell_id$d8d1dc93-44e7-466f-93fe-a9b427368944downstream_cells_mapupstream_cells_map@md_strResourcegetindex$9c31a2c5-e0c7-4891-ad45-7a0e10012bd8precedence_heuristic cell_id$9c31a2c5-e0c7-4891-ad45-7a0e10012bd8downstream_cells_mapupstream_cells_map@md_strgetindex$9643aa92-18fe-4850-a246-8ad134c085e8precedence_heuristic cell_id$9643aa92-18fe-4850-a246-8ad134c085e8downstream_cells_mapupstream_cells_map@md_strgetindex$45a18eea-80c0-4aa4-86d4-c25f90df9303precedence_heuristic cell_id$45a18eea-80c0-4aa4-86d4-c25f90df9303downstream_cells_mapupstream_cells_map@md_strgetindex$85c0b164-b6e0-43f5-a458-939b4fd7dddeprecedence_heuristic cell_id$85c0b164-b6e0-43f5-a458-939b4fd7dddedownstream_cells_mapupstream_cells_map@md_strResourcegetindex$b7244973-218a-4bb6-8a8f-95f84a8fb225precedence_heuristic cell_id$b7244973-218a-4bb6-8a8f-95f84a8fb225downstream_cells_mapupstream_cells_map@md_strResource=>getindex$d058eddf-a8b1-453b-a440-885a3da0a1a0precedence_heuristic cell_id$d058eddf-a8b1-453b-a440-885a3da0a1a0downstream_cells_mapupstream_cells_mapTableOfContents$72412346-af0c-49dc-8c0d-c5b4fe97429bprecedence_heuristic cell_id$72412346-af0c-49dc-8c0d-c5b4fe97429bdownstream_cells_mapupstream_cells_map@md_strResourcegetindex$bf6d18d0-ee99-4a2e-badb-c1f42071903bprecedence_heuristiccell_id$bf6d18d0-ee99-4a2e-badb-c1f42071903bdownstream_cells_mapStaticArraysTestParticlePlutoUIOrdinaryDiffEqupstream_cells_map$05aee60b-fdb4-4856-a1cd-a0dbb470bbb3precedence_heuristic cell_id$05aee60b-fdb4-4856-a1cd-a0dbb470bbb3downstream_cells_mapupstream_cells_map@md_strResourcegetindex$36be1db4-39d5-4067-b42d-53b669bf4ae3precedence_heuristic cell_id$36be1db4-39d5-4067-b42d-53b669bf4ae3downstream_cells_mapupstream_cells_map@md_strgetindex$d50e4350-49ac-4203-870b-903a8801156eprecedence_heuristic cell_id$d50e4350-49ac-4203-870b-903a8801156edownstream_cells_mapupstream_cells_mapBaseBase.Docs.HTML@html_str$45b02635-24ef-4a95-adbd-f6fee4413f21precedence_heuristic cell_id$45b02635-24ef-4a95-adbd-f6fee4413f21downstream_cells_mapupstream_cells_map@md_strResourcegetindex$84dcadfc-4b85-48fc-b7a5-5bbc71d66b7fprecedence_heuristic cell_id$84dcadfc-4b85-48fc-b7a5-5bbc71d66b7fdownstream_cells_mapupstream_cells_map@md_strResourcegetindex$75936fde-d5a7-4920-a4ca-a626a2e90123precedence_heuristic cell_id$75936fde-d5a7-4920-a4ca-a626a2e90123downstream_cells_mapupstream_cells_map@md_strResourcegetindex$0906fe4a-8f2d-4e6d-8726-97bc493da577precedence_heuristic cell_id$0906fe4a-8f2d-4e6d-8726-97bc493da577downstream_cells_mapupstream_cells_map@md_strResourcegetindex$db26a01b-b01c-41dc-b23a-02dbb3867655precedence_heuristic cell_id$db26a01b-b01c-41dc-b23a-02dbb3867655downstream_cells_mapupstream_cells_map@md_strResourcegetindex$47fd7e19-9859-404e-a0e0-be356eef83b0precedence_heuristic cell_id$47fd7e19-9859-404e-a0e0-be356eef83b0downstream_cells_mapupstream_cells_map@md_strResourcegetindex$b3e0e1b5-4f11-41b4-b2e3-fdb36c6c8b20precedence_heuristic cell_id$b3e0e1b5-4f11-41b4-b2e3-fdb36c6c8b20downstream_cells_mapupstream_cells_map@md_strResourcegetindex$96a6b000-b558-11ee-1024-4d2ef81f5ad0precedence_heuristic cell_id$96a6b000-b558-11ee-1024-4d2ef81f5ad0downstream_cells_mapupstream_cells_map@md_strgetindex$021d3bb6-ceb6-460a-a56e-c8b094f3d932precedence_heuristic cell_id$021d3bb6-ceb6-460a-a56e-c8b094f3d932downstream_cells_mapupstream_cells_map@md_strgetindex$e8daa564-73a1-4306-87ab-7334cb44f96dprecedence_heuristic cell_id$e8daa564-73a1-4306-87ab-7334cb44f96ddownstream_cells_mapupstream_cells_map@md_strgetindexcell_execution_order$bf6d18d0-ee99-4a2e-badb-c1f42071903b$d50e4350-49ac-4203-870b-903a8801156e$d058eddf-a8b1-453b-a440-885a3da0a1a0$96a6b000-b558-11ee-1024-4d2ef81f5ad0$75936fde-d5a7-4920-a4ca-a626a2e90123$36be1db4-39d5-4067-b42d-53b669bf4ae3$45a18eea-80c0-4aa4-86d4-c25f90df9303$021d3bb6-ceb6-460a-a56e-c8b094f3d932$e8daa564-73a1-4306-87ab-7334cb44f96d$e73e1a00-3706-42e1-acd1-ebc8b2b95c62$84dcadfc-4b85-48fc-b7a5-5bbc71d66b7f$b7244973-218a-4bb6-8a8f-95f84a8fb225$db26a01b-b01c-41dc-b23a-02dbb3867655$b3e0e1b5-4f11-41b4-b2e3-fdb36c6c8b20$9c31a2c5-e0c7-4891-ad45-7a0e10012bd8$0906fe4a-8f2d-4e6d-8726-97bc493da577$85c0b164-b6e0-43f5-a458-939b4fd7ddde$45b02635-24ef-4a95-adbd-f6fee4413f21$47fd7e19-9859-404e-a0e0-be356eef83b0$d8d1dc93-44e7-466f-93fe-a9b427368944$72412346-af0c-49dc-8c0d-c5b4fe97429b$9643aa92-18fe-4850-a246-8ad134c085e8$081ce252-d083-4f03-a5bc-a4b171d222bc$05aee60b-fdb4-4856-a1cd-a0dbb470bbb3last_hot_reload_timeshortpathtestparticle_202401.jlprocess_statusreadypathJ/home/runner/work/pluto_playground/pluto_playground/testparticle_202401.jlpluto_versionv0.20.13last_save_timeA&~h^cell_order$d50e4350-49ac-4203-870b-903a8801156e$bf6d18d0-ee99-4a2e-badb-c1f42071903b$d058eddf-a8b1-453b-a440-885a3da0a1a0$96a6b000-b558-11ee-1024-4d2ef81f5ad0$75936fde-d5a7-4920-a4ca-a626a2e90123$36be1db4-39d5-4067-b42d-53b669bf4ae3$45a18eea-80c0-4aa4-86d4-c25f90df9303$021d3bb6-ceb6-460a-a56e-c8b094f3d932$e8daa564-73a1-4306-87ab-7334cb44f96d$e73e1a00-3706-42e1-acd1-ebc8b2b95c62$84dcadfc-4b85-48fc-b7a5-5bbc71d66b7f$b7244973-218a-4bb6-8a8f-95f84a8fb225$db26a01b-b01c-41dc-b23a-02dbb3867655$b3e0e1b5-4f11-41b4-b2e3-fdb36c6c8b20$9c31a2c5-e0c7-4891-ad45-7a0e10012bd8$0906fe4a-8f2d-4e6d-8726-97bc493da577$85c0b164-b6e0-43f5-a458-939b4fd7ddde$45b02635-24ef-4a95-adbd-f6fee4413f21$47fd7e19-9859-404e-a0e0-be356eef83b0$d8d1dc93-44e7-466f-93fe-a9b427368944$72412346-af0c-49dc-8c0d-c5b4fe97429b$9643aa92-18fe-4850-a246-8ad134c085e8$081ce252-d083-4f03-a5bc-a4b171d222bc$05aee60b-fdb4-4856-a1cd-a0dbb470bbb3published_objectsnbpkginstall_time_nsTOW뗬instantiatedòinstalled_versionsStaticArrays1.9.1TestParticle0.4.1PlutoUI0.7.55OrdinaryDiffEq6.69.0terminal_outputsStaticArraysX Resolving... === ┌ Warning: Pkg operation failed. Fixing stdlib dependencies and trying again... └ @ GracefulPkg ~/.julia/packages/GracefulPkg/LmTcK/src/apply strategies.jl:92 Resolving... ===  No Changes to `/tmp/jl_T0dgzP/Project.toml`  No Changes to `/tmp/jl_T0dgzP/Manifest.toml` Instantiating... === Precompiling... ===  Activating project at `/tmp/jl_T0dgzP` Precompiling project... 176 dependencies successfully precompiled in 358 seconds. 46 already precompiled. 1 dependency had output during precompilation: ┌ MKL_jll │ Downloading artifact: MKL └ nbpkg_syncX Resolving... === ┌ Warning: Pkg operation failed. Fixing stdlib dependencies and trying again... └ @ GracefulPkg ~/.julia/packages/GracefulPkg/LmTcK/src/apply strategies.jl:92 Resolving... ===  No Changes to `/tmp/jl_T0dgzP/Project.toml`  No Changes to `/tmp/jl_T0dgzP/Manifest.toml` Instantiating... === Precompiling... ===  Activating project at `/tmp/jl_T0dgzP` Precompiling project... 176 dependencies successfully precompiled in 358 seconds. 46 already precompiled. 1 dependency had output during precompilation: ┌ MKL_jll │ Downloading artifact: MKL └ TestParticleX Resolving... === ┌ Warning: Pkg operation failed. Fixing stdlib dependencies and trying again... └ @ GracefulPkg ~/.julia/packages/GracefulPkg/LmTcK/src/apply strategies.jl:92 Resolving... ===  No Changes to `/tmp/jl_T0dgzP/Project.toml`  No Changes to `/tmp/jl_T0dgzP/Manifest.toml` Instantiating... === Precompiling... ===  Activating project at `/tmp/jl_T0dgzP` Precompiling project... 176 dependencies successfully precompiled in 358 seconds. 46 already precompiled. 1 dependency had output during precompilation: ┌ MKL_jll │ Downloading artifact: MKL └ PlutoUIX Resolving... === ┌ Warning: Pkg operation failed. Fixing stdlib dependencies and trying again... └ @ GracefulPkg ~/.julia/packages/GracefulPkg/LmTcK/src/apply strategies.jl:92 Resolving... ===  No Changes to `/tmp/jl_T0dgzP/Project.toml`  No Changes to `/tmp/jl_T0dgzP/Manifest.toml` Instantiating... === Precompiling... ===  Activating project at `/tmp/jl_T0dgzP` Precompiling project... 176 dependencies successfully precompiled in 358 seconds. 46 already precompiled. 1 dependency had output during precompilation: ┌ MKL_jll │ Downloading artifact: MKL └ OrdinaryDiffEqX Resolving... === ┌ Warning: Pkg operation failed. Fixing stdlib dependencies and trying again... └ @ GracefulPkg ~/.julia/packages/GracefulPkg/LmTcK/src/apply strategies.jl:92 Resolving... ===  No Changes to `/tmp/jl_T0dgzP/Project.toml`  No Changes to `/tmp/jl_T0dgzP/Manifest.toml` Instantiating... === Precompiling... ===  Activating project at `/tmp/jl_T0dgzP` Precompiling project... 176 dependencies successfully precompiled in 358 seconds. 46 already precompiled. 1 dependency had output during precompilation: ┌ MKL_jll │ Downloading artifact: MKL └ enabled÷restart_recommended_msgrestart_required_msgbusy_packageswaiting_for_permission,waiting_for_permission_but_probably_disabled«cell_inputs$081ce252-d083-4f03-a5bc-a4b171d222bccell_id$081ce252-d083-4f03-a5bc-a4b171d222bccodemd""" ### Procedures to obtain ion surface flux 1. Obtain magnetic field, density, velocity and current from the unstructured grid. Requires a conversion from BATSRUS unstructured *.out format to VTK *.vtu format. 2. Interpolate MHD quantities onto a Cartesian grid and calculate the derived electric field including the convective term and Hall term. ```math \mathbf{E} = -\sum_\alpha\mathbf{U}_\alpha\times\mathbf{B} + \frac{\mathbf{J}\times\mathbf{B}}{\sum_\alpha q_\alpha n_\alpha} - \frac{\nabla P_e}{\sum_\alpha q_\alpha n_\alpha} + \eta \mathbf{J} ``` This step is performed in ParaView with Python. > It would be better to first convert the basic quantities from spherical coordinates to Cartesian coordinates, then calculate the derived electric field. If we first calculate the electric field and then interpolate to Cartesian coordinates, then around the inner boundary region there might be some artifacts due to interpolations. 3. Read the field info into TestParticle.jl. 4. Tracing with prescribed seeds. ### Potential Issues * In some MHD models (e.g. [Jia+ 2009](https://doi.org/10.1029/2009JA014375)), there is an additional artificial resistivity. * In [Jia+ 2009](https://doi.org/10.1029/2009JA014375)'s MHD model, the energetic population's contribution is included in the pressure term. Whether or not self-consistency is maintained is to be debated. * E field near the inner boundary. """metadatashow_logsèdisabled®skip_as_script«code_folded$e73e1a00-3706-42e1-acd1-ebc8b2b95c62cell_id$e73e1a00-3706-42e1-acd1-ebc8b2b95c62codebegin # Magnetic field B(x) = SA[0.0, 0.0, 1e-8] # [nT] # Electric field E(x) = SA[0.0, 0.0, 0.0] # [V/m] # Initial conditions x0 = [1.0, 0.0, 0.0] # [m] v0 = [0.0, 1.0, 0.1] # [m/s] stateinit = [x0..., v0...] tspan = (0, 18) # [s] # Assemble particle + fields param = prepare(E, B, species=Proton) prob = ODEProblem(trace!, stateinit, tspan, param) # Trace trajectory and save positions & velocities sol = solve(prob, Vern9()) endmetadatashow_logsèdisabled®skip_as_script«code_folded$d8d1dc93-44e7-466f-93fe-a9b427368944cell_id$d8d1dc93-44e7-466f-93fe-a9b427368944codemmd""" ## WIP: Tracing Particles Near Ganymede Ganymede is believed to have both an ionosphere and an exosphere besides the magnetosphere. $(Resource("https://raw.githubusercontent.com/henry2004y/pluto_playground/master/figures/MIE_coupling.png")) Our ultimate goal is to build a self-consistent magnetosphere-ionosphere-exosphere coupling model for Ganymede. """metadatashow_logsèdisabled®skip_as_script«code_folded$9c31a2c5-e0c7-4891-ad45-7a0e10012bd8cell_id$9c31a2c5-e0c7-4891-ad45-7a0e10012bd8codeAmd""" Relativistic B-only tracing: ```math \begin{aligned} \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} &= \frac{q}{\gamma m c}\mathbf{u}\times\mathbf{B} \\ \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} &= \mathbf{u} \end{aligned} ``` Initial particle seeds: random location, $u = c$, pitch angle $\theta = \cos^{-1}(\mu_0),\,\text{where }\mu_0 = 0.1-0.7$, random azimuthal angle. The Lorentz equation in SI units is written as ```math \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t} = \frac{q}{m}\left( \mathbf{v}\times\mathbf{B} + \mathbf{E} \right) ``` It can be normalized to ```math \frac{\mathrm{d}\mathbf{v}^\prime}{\mathrm{d}t^\prime} = \mathbf{v}^\prime\times\mathbf{B}^\prime + \mathbf{E}^\prime ``` with the following transformation ```math \begin{aligned} \mathbf{v} &= \mathbf{v}^\prime V_0 \\ \mathbf{x} &= \mathbf{x}^\prime l_0 \\ t &= t^\prime t_0 \\ \mathbf{B} &= \mathbf{B}^\prime B_0 \\ \mathbf{E} &= \mathbf{E}^\prime E_0 = \mathbf{E}^\prime V_0 B_0 \end{aligned} ``` where $V_0, l_0, t_0, B_0$ are prescribed parameters. After the normalization, we can think of everything in the dimensionless natural units. If the magnetic field is homogeneous and the initial perpendicular velocity is 1, then the gyroradius is 1. MHD solutions are also dimensionless by nature. For the turbulence EM fields, we know their values and the spatial range. To make them more easily understandable, we can convert from the dimensionless MHD to natural units via ```math \begin{aligned} V_0 &= c \\ l_0 &= 4 / nx \\ t_0 &= l_0 / V_0 \\ B_0 &= \mathrm{norm}(B) \end{aligned} ``` and ```julia x = range(xmin, xmax, length=nx) / l₀ y = range(ymin, ymax, length=ny) / l₀ z = range(zmin, zmax, length=nz) / l₀ B ./= B₀ ``` When the X extent is 1, then after the normalization, the x grid resolution becomes 1/4, which is a quarter of the gyroradius. """metadatashow_logsèdisabled®skip_as_script«code_folded$9643aa92-18fe-4850-a246-8ad134c085e8cell_id$9643aa92-18fe-4850-a246-8ad134c085e8codemd""" > The surface flux is very sensitive to the seeding procedures. Placing a plane source in the upstream as in [Plainaki+ 2015, 2020] is deficient in capturing all the possible sources, which ends up in significantly underestimating ion sputtering. ### Stop Conditions * Hitting the surface: $r_\mathrm{final} < 1.07\, R_G$ * Escaping from the outer boundary: $r_\mathrm{final} > 4\,R_G$ """metadatashow_logsèdisabled®skip_as_script«code_folded$45a18eea-80c0-4aa4-86d4-c25f90df9303cell_id$45a18eea-80c0-4aa4-86d4-c25f90df9303codemd""" ## Julia Community Ecosystem - `DifferentialEquations.jl` is one of the earliest and largest packages developed in the community. Among them, OrdinaryDiffEq is a general solver library for ordinary differential equations. - `Interpolations.jl` is the most commonly used high-performance interpolation library in the community. It provides support for different grids, different precisions, and different algorithms. - `StaticArrays.jl` is an optimization library for stack-allocated arrays in the community. It provides memory optimization for small fixed-size arrays that Julia itself does not currently support. - `Meshes.jl` is a new library for different types of grids. - `SpecialFunctions.jl` is a native Julia implementation of special mathematical functions. - `Makie.jl` is a new generation of Julia native graphics library target at GPU rendering. ## Advantages of Julia's Open Source Community - Lego-style combination, avoid reinventing the wheel - Individual abilities are limited - Taking advantage of the vast open-source algorithm library - The function interface at the language level is very convenient to use - Support for generic types - Consideration for performance - Collaboration with strangers - Julia has a large and active community, which makes it easy to find people to collaborate with. This can be a great way to learn new things and get feedback on your work. """metadatashow_logsèdisabled®skip_as_script«code_folded$85c0b164-b6e0-43f5-a458-939b4fd7dddecell_id$85c0b164-b6e0-43f5-a458-939b4fd7dddecodeَmd""" $\mu=0.1,\,l_0=16/nx$ $(Resource("https://raw.githubusercontent.com/henry2004y/pluto_playground/master/figures/CR3_mu0.1_l16.png")) """metadatashow_logsèdisabled®skip_as_script«code_folded$b7244973-218a-4bb6-8a8f-95f84a8fb225cell_id$b7244973-218a-4bb6-8a8f-95f84a8fb225code$md""" ### Magnetospheric Ion Motion - Spitzer orbit - Gradient and curvature drifts - Orbits near the tail reconnection site $(Resource("https://raw.githubusercontent.com/henry2004y/TestParticlePresentation/master/electron_fixedtime_pm0.5RE_EGI_t1298.mp4", :autoplay => "", :loop => "")) """metadatashow_logsèdisabled®skip_as_script«code_folded$d058eddf-a8b1-453b-a440-885a3da0a1a0cell_id$d058eddf-a8b1-453b-a440-885a3da0a1a0codeTableOfContents()metadatashow_logsèdisabled®skip_as_script«code_folded$72412346-af0c-49dc-8c0d-c5b4fe97429bcell_id$72412346-af0c-49dc-8c0d-c5b4fe97429bcodemd""" ### Our Model We use ideal MHD G2 spherical output as a starting point. - G2 ideal MHD - G2 Hall MHD - G8 ideal MHD - G8 Hall MHD ### Seeding An important question is how to set seeds. * Thermal population: using MHD values * Energetic population: using observed differential flux, 10 logarithmic bins for sampling $(Resource("https://agupubs.onlinelibrary.wiley.com/cms/asset/818da498-ee6b-4e80-9136-46d032b32092/jgra54295-fig-0002-m.jpg")) """metadatashow_logsèdisabled®skip_as_script«code_folded$bf6d18d0-ee99-4a2e-badb-c1f42071903bcell_id$bf6d18d0-ee99-4a2e-badb-c1f42071903bcodebbegin using TestParticle using OrdinaryDiffEq using StaticArrays using PlutoUI endmetadatashow_logsèdisabled®skip_as_script«code_folded$05aee60b-fdb4-4856-a1cd-a0dbb470bbb3cell_id$05aee60b-fdb4-4856-a1cd-a0dbb470bbb3code]md""" ### Preliminary Results * Thermal O+ population $(Resource("https://raw.githubusercontent.com/henry2004y/pluto_playground/master/figures/impact_flux_G2_thermal_O_idealMHD.png")) * Energetic H+ population $(Resource("https://raw.githubusercontent.com/henry2004y/pluto_playground/master/figures/impact_flux_G2_idealMHD_H_bin1-5-9.png")) """metadatashow_logsèdisabled®skip_as_script«code_folded$36be1db4-39d5-4067-b42d-53b669bf4ae3cell_id$36be1db4-39d5-4067-b42d-53b669bf4ae3codemd""" ## Motion of Charged Particles ### Non-relativistic case ```math m\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t} = q(\mathbf{E} + \mathbf{v}\times\mathbf{B}) \tag{1} ``` ### Relativistic case The four-dimensional momentum equation satisfying Lorentz covariance is ```math \frac{\mathrm{d}p^\mu}{\mathrm{d}\tau} = \gamma qF^{\mu \nu}v_\nu \tag{2} ``` where $\gamma = 1/\sqrt{1-v^2/c^2}$ is the Lorentz factor， $\tau$ is the proper time， $p^\mu = \gamma m v^\mu$ is the four-dimensional momentum, and $F^{\mu\nu}$ is the electromagnetic tensor。 There are two forms that can be used in numerical computation after simplification: ```math \frac{\mathrm{d}(\gamma v^i)}{\mathrm{d}t} = \frac{q}{m}(E^i+\epsilon^{ijk}v_j B_k) \tag{3} ``` or ```math \frac{\mathrm{d}v^i}{\mathrm{d}t} = \frac{q}{m \gamma^3}(E^i+\epsilon^{ijk}v_j B_k) \tag{4} ``` Currently we only include Eq.(4) in TestParticle.jl. This needs further testing. """metadatashow_logsèdisabled®skip_as_script«code_folded$d50e4350-49ac-4203-870b-903a8801156ecell_id$d50e4350-49ac-4203-870b-903a8801156ecode2html""metadatashow_logsèdisabled®skip_as_script«code_folded$45b02635-24ef-4a95-adbd-f6fee4413f21cell_id$45b02635-24ef-4a95-adbd-f6fee4413f21codeٌmd""" $\mu=0.7,\,l_0=4/nx$ $(Resource("https://raw.githubusercontent.com/henry2004y/pluto_playground/master/figures/CR1_mu0.7_l4.png")) """metadatashow_logsèdisabled®skip_as_script«code_folded$84dcadfc-4b85-48fc-b7a5-5bbc71d66b7fcell_id$84dcadfc-4b85-48fc-b7a5-5bbc71d66b7fcodemd""" ## Demos ### Electron and Ion Motion $(Resource("https://raw.githubusercontent.com/henry2004y/TestParticlePresentation/master/double_trajectories_static_B.png")) ### Grad-B and Curvature Drift $(Resource("https://henry2004y.github.io/TestParticle.jl/dev/examples/basics/demo_curvature_B-4.png")) ### Finite-Larmor-Radius Effect $(Resource("https://henry2004y.github.io/TestParticle.jl/dev/examples/basics/demo_FLR-4.png")) ### Polarization Drift $(Resource("https://henry2004y.github.io/TestParticle.jl/dev/examples/basics/demo_polarization_drift-4.png")) ### Adiabatic motion $(Resource("https://raw.githubusercontent.com/henry2004y/TestParticlePresentation/master/adiabatic_trajectory_dipole.png")) """metadatashow_logsèdisabled®skip_as_script«code_folded$75936fde-d5a7-4920-a4ca-a626a2e90123cell_id$75936fde-d5a7-4920-a4ca-a626a2e90123code md""" ## What Are Test Particles? A test particle model describes non-self-consistent field-particle interactions. A complete description of the particle phase space trajectory not only needs to reflect the effect of the field on the particle, but also needs to reflect the effect of the particle on the field. **Field** <--> **Particle** However, as an approximate treatment, we can ignore the feedback of a small number of particles on the field and only focus on the effect of the field on the particle. **Field** --> **Particle** $(Resource("https://www.schoolphysics.co.uk/age14-16/Electricity%20and%20magnetism/Current%20electricity/text/Particles_and_fields/images/2.png")) ## Applications of Test Particle Models - Studying charged particle trajectories under EM fields. - Explaining collective behaviors of particles controlled by fields. - Substituting expensive laboratory experiments. $(Resource("https://raw.githubusercontent.com/henry2004y/TestParticlePresentation/master/single_trajectory_static_B.png")) """metadatashow_logsèdisabled®skip_as_script«code_folded$0906fe4a-8f2d-4e6d-8726-97bc493da577cell_id$0906fe4a-8f2d-4e6d-8726-97bc493da577codemd""" ### Scales in the System * Numerical dissipation scale * Injection scale $\mu=0.1,\,l_0=4/nx$ $(Resource("https://raw.githubusercontent.com/henry2004y/pluto_playground/master/figures/CR3_mu0.1_l4.png")) """metadatashow_logsèdisabled®skip_as_script«code_folded$db26a01b-b01c-41dc-b23a-02dbb3867655cell_id$db26a01b-b01c-41dc-b23a-02dbb3867655codebmd""" ## Demo: Tracing Electrons in Mars' Crustal Field $(Resource("https://raw.githubusercontent.com/henry2004y/pluto_playground/master/figures/mars_electron_dispersion.png")) Goal: Given Mars' crustal field and the locations of the observed electrons, can we determine their source location? Unknown factors: * Exact pitch angles * Electric field $(Resource("https://raw.githubusercontent.com/henry2004y/pluto_playground/master/figures/mars_crustal_field_pass_count.png")) $(Resource("https://raw.githubusercontent.com/henry2004y/pluto_playground/master/figures/electron_trajectory_3eV_10eV.png")) """metadatashow_logsèdisabled®skip_as_script«code_folded$47fd7e19-9859-404e-a0e0-be356eef83b0cell_id$47fd7e19-9859-404e-a0e0-be356eef83b0codeَmd""" $\mu=0.7,\,l_0=16/nx$ $(Resource("https://raw.githubusercontent.com/henry2004y/pluto_playground/master/figures/CR1_mu0.7_l16.png")) """metadatashow_logsèdisabled®skip_as_script«code_folded$b3e0e1b5-4f11-41b4-b2e3-fdb36c6c8b20cell_id$b3e0e1b5-4f11-41b4-b2e3-fdb36c6c8b20code(md""" ## Demo: Tracing Particles in MHD Turbulence $(Resource("https://raw.githubusercontent.com/henry2004y/pluto_playground/master/figures/turbulence_mhd_bx_cut.png")) ```math \begin{aligned} B_\mathrm{in-plane} &= 1 \\ B_z &= 4 \\ U_\mathrm{in-plane} &= 1 \\ U_z &= 0 \\ \end{aligned} ``` """metadatashow_logsèdisabled®skip_as_script«code_folded$96a6b000-b558-11ee-1024-4d2ef81f5ad0cell_id$96a6b000-b558-11ee-1024-4d2ef81f5ad0codeTmd""" # TestParticle.jl: A New Tool for An Old Problem 2024/01/31 Hongyang Zhou """metadatashow_logsèdisabled®skip_as_script«code_folded$021d3bb6-ceb6-460a-a56e-c8b094f3d932cell_id$021d3bb6-ceb6-460a-a56e-c8b094f3d932code(md""" ## The Goal of TestParticle.jl 1. Simplicity: Users only need to provide field information and particle initial conditions 2. Generality: Compatible with general models, can run independently 3. Efficiency: The efficiency of C 4. Visuality: Connected to powerful visualization tools """metadatashow_logsèdisabled®skip_as_script«code_folded$e8daa564-73a1-4306-87ab-7334cb44f96dcell_id$e8daa564-73a1-4306-87ab-7334cb44f96dcode md""" ## Design ### DifferentialEquations.jl Interface ```julia param = prepare(E, B, species) prob = ODEProblem(trace!, stateinit, tspan, param) sol = solve(prob, algorithm) ``` - Analytical or numerical field? - Contains external field (e.g. gravity)? - Normalized unit system? - Time-dependent? - Relativistic? By taking advantage of Julia's multiple dispatch property, we come up with 4 ODE functions: - `trace!` - `trace_relativistic!` - `trace_normalized!` - `trace_relativistic_normalized!` ### Energy Conservation We need to make sure that the energy is conserved in a B-only tracing for validity. In classical mechanics, Hamiltonian's equations naturally have that the solutions reside on a symplectic manifold in phase space, with the natural splitting being the position and momentum components. _For the true Hamiltonian solution, that phase space path is constant energy_. For example, the most commonly used solver in DifferentialEquations.jl, `Tsit5()` (adaptive 4/5th order Runge-Kutta method), does not conserve energy well. The following shows the test result of a proton in a static magnetic field gyrating for 304 cycles. All adaptive timestepping schemes from DifferentialEquations.jl use a relative tolerence of 1e-3 and an absolute tolerence of 1e-6. | Scheme | Energy difference ratio | Time [ms] | Memory [MiB] | fixed time step | |-----------|-------------|---------|-----------|-----| | **Boris** | 0% | 0.2 | 0.047 | true | | **ImplicitMidpoint** | 0% | 6.1 | 2.974 | true | | **SSPSDIRK2** | 0% | 11.7 | 4.090 | true | | Vern8 | 0.17% | 0.9 | 1.707 | false | | RosenbrockW6S4OS | 0.23% | 11.5 | 3.464 | true | | Vern9 | -0.59% | 1.1 | 1.662 | false | | Trapezoid | 1.14% | 20.2 | 5.088 | false | | Vern6 | -5.36% | 1.1 | 2.036 | false | | Tsit5 | 20.59% | 1.2 | 2.309 | false | | BS3 | -89.73% | 2.6 | 3.503 | false | ### Native Boris Solver - `trace_trajectory` The classical Boris method is not a general ODE solver, but designed specifically for the Lorentz equation with the symmetry of the EM fields. ```julia param = prepare(E, B, species=Electron) prob = TraceProblem(stateinit, tspan, param) sol = trace_trajectory(prob; dt, savestepinterval) ``` This method is 2nd order accurate in space. A fixed `dt` is required to conserve phase space volume, even though it is not a symplectic method ([Qin+ 2013](https://doi.org/10.1063/1.4818428)). It is the de facto algorithm for particle pusher after being proposed in the 1960s. However, the phase error may still be an issue. ### Implementation: Analytic or Numerical Fields? - Analytic --> `Field(F)` - Numerical --> Needs interpolation ```julia function getinterp(A, gridx, gridy, gridz, order::Int=1, bc::Int=1) @assert size(A,1) == 3 && ndims(A) == 4 "Inconsistent 3D force field and grid!" Ax = @view A[1,:,:,:] Ay = @view A[2,:,:,:] Az = @view A[3,:,:,:] itpx, itpy, itpz = _getinterp(Ax, Ay, Az, order, bc) interpx = scale(itpx, gridx, gridy, gridz) interpy = scale(itpy, gridx, gridy, gridz) interpz = scale(itpz, gridx, gridy, gridz) # Return field value at a given location. function get_field(xu) r = @view xu[1:3] return SA[interpx(r...), interpy(r...), interpz(r...)] end return get_field end ``` """metadatashow_logsèdisabled®skip_as_script«code_foldedënotebook_id$d09ca6a8-76bc-11f0-1aa6-e505a09f9b3fin_temp_dir¨metadata