henry2004y / TestParticle.jl

Test particle tracing in electromagnetic field
https://henry2004y.github.io/TestParticle.jl/dev/
MIT License
16 stars 5 forks source link

Save the final timepoint for Boris method #194

Open Beforerr opened 1 month ago

Beforerr commented 1 month ago

Right now, it is not possible to save the last step if you don't want to take any intermediate steps. It would be nice to have a kwarg like save_end: Denotes whether the final timepoint is forced to be saved, regardless of the other saving settings.

See https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts

henry2004y commented 1 month ago

Thanks for the feature request! I guess this should be pretty easy to add in https://github.com/henry2004y/TestParticle.jl/blob/master/src/pusher.jl

henry2004y commented 1 month ago

I found the following logic in OrdinaryDiffEq.jl regarding save_end and save_start:

https://github.com/SciML/OrdinaryDiffEq.jl/blob/873f4400e8242e0f545adf5f240f04d507ae0d11/lib/OrdinaryDiffEqCore/src/solve.jl#L285

henry2004y commented 1 month ago

Hi @Beforerr, do you have a minimal working example for me to test on?

Beforerr commented 1 month ago

Something like that? I modified from your example with save_everystep and save_end keywords ( it is not working)

# ---
# title: Boris method save keywords
# id: demo_boris_save
# ---

using TestParticle
using StaticArrays
using OrdinaryDiffEq

const Bmag = 0.01
uniform_B(x) = SA[0.0, 0.0, Bmag]
uniform_E(x) = SA[0.0, 0.0, 0.0]

## Reference parameters
const tperiod = 2π / (abs(param[1]) * sqrt(sum(x -> x^2, param[3]([0.0,0.0,0.0], 0.0))))
const rL = sqrt(v0[1]^2 + v0[2]^2 + v0[3]^2) / (abs(param[1]) * Bmag)
const invrL = 1 / rL;

tspan = (0.0, 10 * tperiod)
dt = tperiod / 4

x0 = [0.0, 0.0, 0.0]
v0 = [0.0, 1e5, 0.0]
u0 = [x0..., v0...]
u0s = [u0 for _ in 1:1000]

param = prepare(uniform_E, uniform_B, species=Electron)
prob_func = (prob, i, repeat = nothing) -> remake(prob, u0 = u0s[i])
prob = TraceProblem(u0, tspan, param; prob_func)
sol_boris = TestParticle.solve(prob; trajectories=length(u0s), dt, save_everystep = false, save_end = true)