PrincetonUniversity / SPEC

The Stepped-Pressure Equilibrium Code, an advanced MRxMHD equilibrium solver.
https://princetonuniversity.github.io/SPEC/
GNU General Public License v3.0
24 stars 4 forks source link

Speed up field line tracing #79

Open krystophny opened 5 years ago

krystophny commented 5 years ago

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 commented 5 years ago

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 notifications@github.com 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 commented 5 years ago

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 commented 5 years ago

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 < notifications@github.com> 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 commented 5 years ago

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 commented 5 years ago

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.