Speed up field line tracing

[krystophny]

After discussion at ISHW with @jonathanschilling we had the idea to speed up field line tracing by using symplectic integration.

Another, possibly even faster and simpler alternative would be to construct an interpolation of the single-turn map between Poincaré cuts and then evaluate the map for the number of turns. Construction requires 1e4-1e6 field lines to be traced over a single toroidal turn with the existing Runge-Kutta scheme and can be done in parallel.

See e.g. Kasilov et al, Phys. Plasmas 9, 3508 (2002); https://doi.org/10.1063/1.1493793 which is already much more complicated as particle orbits are mapped instead of field lines. This is also related to #47 as on can start tracing at arbitrary theta based on such a map.

[acerfon]

Dear all, I do not think symplectic integration would lead to faster field line tracing. For a given accuracy, a non-symplectic ODE solver will be faster than a symplectic one - I believe it is almost a theorem. If you are guaranteed to have flux surfaces, such as in tokamaks, then symplectic integration makes a lot of sense, because you can afford to lose a lot of accuracy while still being sure to remain on the flux surface thanks to the symplectic nature of the solver. However, for general stellarator magnetic fields, it is relatively easy to be much more accurate with a non-symplectic solver than with a symplectic one. Have you thought of spectral deferred correction (SDC) schemes? I believe they are close to being the state-of-the-art ODE solvers for these tasks; this is what we use for field line tracing at Courant. Regards, Antoine

…

On Thu, Oct 10, 2019 at 3:33 AM Christopher Albert ***@***.***> wrote: After discussion at ISHW with @jonathanschilling <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_jonathanschilling&d=DwMFaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=meEavcbFdc0ufv1PWzptlA&m=U4XSQqsKcZFfE4ZqiBIVAJrxW3jlnGk5aCNFm1lk8_I&s=V5pegu5wevBx7AfAqQ75y9MfFwFL5np7hMYslR_Mcq8&e=> we had the idea to speed up field line tracing by using symplectic integration. Another, possibly even faster and simpler alternative would be to construct an interpolation of the single-turn map between Poincaré cuts and then evaluate the map for the number of turns. Construction requires 1e4-1e6 field lines to be traced over a single toroidal turn with the existing Runge-Kutta scheme and can be done in parallel. See e.g. Kasilov et al, Phys. Plasmas 9, 3508 (2002); https://doi.org/10.1063/1.1493793 <https://urldefense.proofpoint.com/v2/url?u=https-3A__doi.org_10.1063_1.1493793&d=DwMFaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=meEavcbFdc0ufv1PWzptlA&m=U4XSQqsKcZFfE4ZqiBIVAJrxW3jlnGk5aCNFm1lk8_I&s=4w_3Lr0du9IWvgmdLxvEL19hI6RgU3qt8wZvR7eOhus&e=> which is already much more complicated as particle orbits are mapped instead of field lines. This is also related to #47 <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_PrincetonUniversity_SPEC_issues_47&d=DwMFaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=meEavcbFdc0ufv1PWzptlA&m=U4XSQqsKcZFfE4ZqiBIVAJrxW3jlnGk5aCNFm1lk8_I&s=cp4bl3TuENHq9pAI6RfOIDAcJ8QdrbKgNW1j8N5xc5c&e=> as on can start tracing at arbitrary theta based on such a map. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_PrincetonUniversity_SPEC_issues_79-3Femail-5Fsource-3Dnotifications-26email-5Ftoken-3DAB5QJUOP4OMNXVUSUCWFQOLQN3LCVA5CNFSM4I7I77BKYY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4HQ3CKZQ&d=DwMFaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=meEavcbFdc0ufv1PWzptlA&m=U4XSQqsKcZFfE4ZqiBIVAJrxW3jlnGk5aCNFm1lk8_I&s=zdfBC-7S-fxtuNUPwSpmnV9FRSQqZ4kxqkOa6R2OUTg&e=>, or unsubscribe <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AB5QJUMGFHMC4KXF2QVFNXTQN3LCVANCNFSM4I7I77BA&d=DwMFaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=meEavcbFdc0ufv1PWzptlA&m=U4XSQqsKcZFfE4ZqiBIVAJrxW3jlnGk5aCNFm1lk8_I&s=bgreTEaTf3piqaq-hJ8ER8JSZdZ3_VeGbHaE3Gu21Gk&e=> .

[krystophny]

Dear Antoine,

thanks for your fast reply. Whether symplectic or non-symplectic schemes are faster at a given accuracy depends mostly on how long you trace the orbit. For alpha particles in stellarators that are traced over 100.000s to millions of toroidal turns we saw that low-order symplectic schemes outperform RK4/5 by factor 3-5 in computation time. However, if you need only 10.000 turns it is likely that adaptive schemes are faster. Anyway, I would propose to use direct integration only for a single turn and then immediately switch to an interpolated mapping, which will be orders of magnitude faster to evaluate than any numerical integration. We have a working Fortran code for that in Graz and I will play a bit with it and let you know if it is as useful as I think it is.

Best,

Chris

[acerfon]

Hi Chris, I fully agree with your single turn/interpolation idea. I also agree that RK4/5 is pretty bad. But then again, do not discard Spectral Deferred Correction. The best one can do, in my opinion, and Vladimir Rokhlin's, and Jon Wilkening's, and a few other not too shabby numerical analysts (which I am not). Regards, Antoine

…

On Thu, Oct 10, 2019 at 10:39 AM Christopher Albert < ***@***.***> wrote: Dear Antoine, thanks for your fast reply. Whether symplectic or non-symplectic schemes are faster at a given accuracy depends mostly on how long you trace the orbit. For alpha particles in stellarators that are traced over 100.000s to millions of toroidal turns we saw that low-order symplectic schemes outperform RK4/5 by factor 3-5 in computation time. However, if you need only 10.000 turns it is likely that adaptive schemes are faster. Anyway, I would propose to use direct integration only for a single turn and then immediately switch to an interpolated mapping, which will be orders of magnitude faster to evaluate than any numerical integration. We have a working Fortran code for that in Graz and I will play a bit with it and let you know if it is as useful as I think it is. Best, Chris — You are receiving this because you commented. Reply to this email directly, view it on GitHub <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_PrincetonUniversity_SPEC_issues_79-3Femail-5Fsource-3Dnotifications-26email-5Ftoken-3DAB5QJUJYVQG2DJQWTMF2IVDQN45BVA5CNFSM4I7I77BKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEA4TCVY-23issuecomment-2D540619095&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=meEavcbFdc0ufv1PWzptlA&m=FyNbZFmVM9_GerrAN45QQ4bx8NycxYrkO2fWbO93uvQ&s=PgwTV774ZL6VO_q8CTt3Bj5GalUus6DiK3z5KR5JCM8&e=>, or unsubscribe <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AB5QJUMPFFRC7DTGEBTH2STQN45BVANCNFSM4I7I77BA&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=meEavcbFdc0ufv1PWzptlA&m=FyNbZFmVM9_GerrAN45QQ4bx8NycxYrkO2fWbO93uvQ&s=PJ6O7fi2ecJdE7lzHCFwsY3B2uxcVHCPLyhUVN0qIFs&e=> .

[krystophny]

Hi Antoine,

sure - didn't mean to dismiss the SDC methods and I would also call myself rather a user than an expert with regard to numerical integrators.

Best,
Chris

[jloizu]

This is also related to #47 as on can start tracing at arbitrary theta based on such a map.
I just want to say that @abaillod has already implemented a small change in the poincare routines so that one can choose in the input file the starting poloidal angle for the field-line tracing. It is currently only in his branch #75 but should be merged into master soon.