Integrators.H

Go to the documentation of this file.

/* Copyright 2022-2023 The Regents of the University of California, through Lawrence
 *           Berkeley National Laboratory (subject to receipt of any required
 *           approvals from the U.S. Dept. of Energy). All rights reserved.
 *
 * This file is part of ImpactX.
 *
 * Authors: Chad Mitchell, Axel Huebl
 * License: BSD-3-Clause-LBNL
 */
#ifndef IMPACTX_INTEGRATORS_H_
#define IMPACTX_INTEGRATORS_H_
 
#include <AMReX_Extension.H>  // for AMREX_RESTRICT
#include <AMReX_REAL.H>       // for ParticleReal
 
 
namespace impactx::integrators
{

    template <typename T_Element>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void symp2_integrate (
        RefPart & refpart,
        amrex::ParticleReal const zin,
        amrex::ParticleReal const zout,
        int const nsteps,
        T_Element const & element
    )
    {
        using namespace amrex::literals; // for _rt and _prt
 
        // initialize numerical integration parameters
        amrex::ParticleReal const dz = (zout - zin) / static_cast<amrex::ParticleReal>(nsteps);
        amrex::ParticleReal const tau1 = dz * 0.5_prt;
        amrex::ParticleReal const tau2 = dz;
 
        // initialize the value of the independent variable
        amrex::ParticleReal zeval = zin;
 
        // loop over integration steps
        for(int j=0; j < nsteps; ++j)
        {
            element.map1(tau1,refpart,zeval);
            element.map2(tau2,refpart,zeval);
            element.map1(tau1,refpart,zeval);
        }
    }

    template <typename T_Element>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void symp2_integrate_split3 (
        RefPart & refpart,
        amrex::ParticleReal const zin,
        amrex::ParticleReal const zout,
        int const nsteps,
        T_Element const & element
    )
    {
        using namespace amrex::literals; // for _rt and _prt
 
        // initialize numerical integration parameters
        amrex::ParticleReal const dz = (zout - zin) / static_cast<amrex::ParticleReal>(nsteps);
        amrex::ParticleReal const tau1 = dz * 0.5_prt;
        amrex::ParticleReal const tau2 = dz * 0.5_prt;
        amrex::ParticleReal const tau3 = dz;
 
        // initialize the value of the independent variable
        amrex::ParticleReal zeval = zin;
 
        // loop over integration steps
        for(int j=0; j < nsteps; ++j)
        {
            element.map1(tau1,refpart,zeval);
            element.map2(tau2,refpart,zeval);
            element.map3(tau3,refpart,zeval);
            element.map2(tau2,refpart,zeval);
            element.map1(tau1,refpart,zeval);
        }
    }

    template <typename T_Element>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void symp4_integrate (
        RefPart & refpart,
        amrex::ParticleReal const zin,
        amrex::ParticleReal const zout,
        int const nsteps,
        T_Element const & element
    )
    {
        using namespace amrex::literals; // for _rt and _prt
 
        // initialize numerical integration parameters
        amrex::ParticleReal const dz = (zout - zin) / static_cast<amrex::ParticleReal>(nsteps);
        constexpr amrex::ParticleReal alpha = -0.2599210498948732;  // C++26: 1_prt - std::pow(2_prt,1_prt/3_prt);
        amrex::ParticleReal const tau2 = dz / (1_prt + alpha);
        amrex::ParticleReal const tau1 = tau2*0.5_prt;
        amrex::ParticleReal const tau3 = alpha*tau1;
        amrex::ParticleReal const tau4 = (alpha - 1_prt)*tau2;
 
        // initialize the value of the independent variable
        amrex::ParticleReal zeval = zin;
 
        // loop over integration steps
        for (int j=0; j < nsteps; ++j)
        {
            element.map1(tau1,refpart,zeval);
            element.map2(tau2,refpart,zeval);
            element.map1(tau3,refpart,zeval);
            element.map2(tau4,refpart,zeval);
            element.map1(tau3,refpart,zeval);
            element.map2(tau2,refpart,zeval);
            element.map1(tau1,refpart,zeval);
        }
    }

    template <typename T_Real, typename T_Element>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void symp2_integrate_particle (
        amrex::SmallVector<T_Real, 6, 1> & particle,
        amrex::ParticleReal const zin,
        amrex::ParticleReal const zout,
        int const nsteps,
        T_Element const & element
    )
    {
        using namespace amrex::literals; // for _rt and _prt
 
        // initialize numerical integration parameters
        amrex::ParticleReal const dz = (zout - zin) / static_cast<amrex::ParticleReal>(nsteps);
        amrex::ParticleReal const tau1 = dz * 0.5_prt;
        amrex::ParticleReal const tau2 = dz;
 
        // initialize the value of the independent variable
        amrex::ParticleReal zeval = zin;
 
        // loop over integration steps
        for(int j=0; j < nsteps; ++j)
        {
            element.map1(tau1,particle,zeval);
            element.map2(tau2,particle,zeval);
            element.map1(tau1,particle,zeval);
        }
    }
 

    template <typename T_Real, typename T_Element>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void symp4_integrate_particle (
        amrex::SmallVector<T_Real, 6, 1> & particle,
        amrex::ParticleReal const zin,
        amrex::ParticleReal const zout,
        int const nsteps,
        T_Element const & element
    )
    {
        using namespace amrex::literals; // for _rt and _prt
 
        // initialize numerical integration parameters
        amrex::ParticleReal const dz = (zout - zin) / static_cast<amrex::ParticleReal>(nsteps);
        constexpr amrex::ParticleReal alpha = -0.2599210498948732;  // C++26: 1_prt - std::pow(2_prt,1_prt/3_prt);
        amrex::ParticleReal const tau2 = dz/(1_prt + alpha);
        amrex::ParticleReal const tau1 = tau2*0.5_prt;
        amrex::ParticleReal const tau3 = alpha*tau1;
        amrex::ParticleReal const tau4 = (alpha - 1_prt)*tau2;
 
        // initialize the value of the independent variable
        amrex::ParticleReal zeval = zin;
 
        // loop over integration steps
        for (int j=0; j < nsteps; ++j)
        {
            element.map1(tau1,particle,zeval);
            element.map2(tau2,particle,zeval);
            element.map1(tau3,particle,zeval);
            element.map2(tau4,particle,zeval);
            element.map1(tau3,particle,zeval);
            element.map2(tau2,particle,zeval);
            element.map1(tau1,particle,zeval);
        }
    }
 

    template <typename T_Real, typename T_Element>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void symp6_integrate_particle (
        amrex::SmallVector<T_Real, 6, 1> & particle,
        amrex::ParticleReal const zin,
        amrex::ParticleReal const zout,
        int const nsteps,
        T_Element const & element
    )
    {
        using namespace amrex::literals; // for _rt and _prt
 
        // initialize numerical integration parameters
        amrex::ParticleReal const dz = (zout - zin) / static_cast<amrex::ParticleReal>(nsteps);
        amrex::ParticleReal const w0 = 1.31518632068391121888424973;
        amrex::ParticleReal const w1 = -1.17767998417887100694641568;
        amrex::ParticleReal const w2 = 0.23557321335935813368479318;
        amrex::ParticleReal const w3 = 0.78451361047755726381949763;
 
        amrex::ParticleReal const tau1 = w3 * dz / 2_prt;
        amrex::ParticleReal const tau2 = w3 * dz;
        amrex::ParticleReal const tau3 = (w3 + w2) * dz / 2_prt;
        amrex::ParticleReal const tau4 = w2 * dz;
        amrex::ParticleReal const tau5 = (w2 + w1) * dz / 2_prt;
        amrex::ParticleReal const tau6 = w1 * dz;
        amrex::ParticleReal const tau7 = (w1 + w0) * dz / 2_prt;
        amrex::ParticleReal const tau8 = w0 * dz;
 
        // initialize the value of the independent variable
        amrex::ParticleReal zeval = zin;
 
        // loop over integration steps
        for (int j=0; j < nsteps; ++j)
        {
            element.map1(tau1,particle,zeval);
            element.map2(tau2,particle,zeval);
            element.map1(tau3,particle,zeval);
            element.map2(tau4,particle,zeval);
            element.map1(tau5,particle,zeval);
            element.map2(tau6,particle,zeval);
            element.map1(tau7,particle,zeval);
            element.map2(tau8,particle,zeval);
            element.map1(tau7,particle,zeval);
            element.map2(tau6,particle,zeval);
            element.map1(tau5,particle,zeval);
            element.map2(tau4,particle,zeval);
            element.map1(tau3,particle,zeval);
            element.map2(tau2,particle,zeval);
            element.map1(tau1,particle,zeval);
        }
    }

    template <typename T_Real, typename T_Element>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void symp2_integrate_particle_spin_step (
        amrex::SmallVector<T_Real, 6, 1> & particle,
        amrex::SmallVector<T_Real, 3, 1> & particle_spin,
        amrex::ParticleReal const tau,
        amrex::ParticleReal & zeval,
        T_Element const & element
    )
    {
        using namespace amrex::literals; // for _rt and _prt
 
        // initialize numerical integration parameters
        amrex::ParticleReal const tau1 = tau*0.5_prt;
        amrex::ParticleReal const tau2 = tau*0.5_prt;
        amrex::ParticleReal const tau3 = tau;
 
        // a single second-order step
        element.map1(tau1,particle,zeval);
        element.map2(tau2,particle,zeval);
        element.map3(tau3,particle,particle_spin,zeval);
        element.map2(tau2,particle,zeval);
        element.map1(tau1,particle,zeval);
    }

    template <typename T_Real, typename T_Element>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void symp4_integrate_particle_spin_step (
        amrex::SmallVector<T_Real, 6, 1> & particle,
        amrex::SmallVector<T_Real, 3, 1> & particle_spin,
        amrex::ParticleReal const tau,
        amrex::ParticleReal & zeval,
        T_Element const & element
    )
    {
        using namespace amrex::literals; // for _rt and _prt
 
        // initialize numerical integration parameters (Yoshida)
        constexpr amrex::ParticleReal a = 1.2599210498948732;  // C++26 std::pow(2_prt,1_prt/3_prt);
        constexpr amrex::ParticleReal z0 = 1_prt / (2_prt - a);
        constexpr amrex::ParticleReal z1 = -a / (2_prt - a);
        amrex::ParticleReal const tau1 = tau*z0;
        amrex::ParticleReal const tau2 = tau*z1;
 
        // a single fourth-order step
        symp2_integrate_particle_spin_step(particle,particle_spin,tau1,zeval,element);
        symp2_integrate_particle_spin_step(particle,particle_spin,tau2,zeval,element);
        symp2_integrate_particle_spin_step(particle,particle_spin,tau1,zeval,element);
    }
 

    template <typename T_Real, typename T_Element>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void symp6_integrate_particle_spin_step (
        amrex::SmallVector<T_Real, 6, 1> & particle,
        amrex::SmallVector<T_Real, 3, 1> & particle_spin,
        amrex::ParticleReal const tau,
        amrex::ParticleReal & zeval,
        T_Element const & element
    )
    {
        using namespace amrex::literals; // for _rt and _prt
 
        // initialize numerical integration parameters (Yoshida)
        constexpr amrex::ParticleReal a = 1.148698354997035;  // C++26 std::pow(2_prt,1_prt/5_prt);
        constexpr amrex::ParticleReal z0 = 1_prt / (2_prt - a);
        constexpr amrex::ParticleReal z1 = -a / (2_prt - a);
        amrex::ParticleReal const tau1 = tau * z0;
        amrex::ParticleReal const tau2 = tau * z1;
 
        // a single sixth-order step
        symp4_integrate_particle_spin_step(particle,particle_spin,tau1,zeval,element);
        symp4_integrate_particle_spin_step(particle,particle_spin,tau2,zeval,element);
        symp4_integrate_particle_spin_step(particle,particle_spin,tau1,zeval,element);
    }
 

    template <typename T_Real, typename T_Element>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void symp_integrate_particle_spin (
        amrex::SmallVector<T_Real, 6, 1> & particle,
        amrex::SmallVector<T_Real, 3, 1> & particle_spin,
        amrex::ParticleReal const zin,
        amrex::ParticleReal const zout,
        int const nsteps,
        int int_order,
        T_Element const & element
    )
    {
        using namespace amrex::literals; // for _rt and _prt
 
        // initialize numerical integration parameters
        amrex::ParticleReal const dz = (zout - zin) / static_cast<amrex::ParticleReal>(nsteps);
 
        // initialize the value of the independent variable
        amrex::ParticleReal zeval = zin;
 
        // loop over integration steps
        for(int j=0; j < nsteps; ++j)
        {
            if (int_order == 2) {
                symp2_integrate_particle_spin_step(particle,particle_spin,dz,zeval,element);
            } else if (int_order == 4) {
                symp4_integrate_particle_spin_step(particle,particle_spin,dz,zeval,element);
            } else if (int_order == 6) {
                symp6_integrate_particle_spin_step(particle,particle_spin,dz,zeval,element);
            } else {
                symp2_integrate_particle_spin_step(particle,particle_spin,dz,zeval,element);
            }
        }
    }
 
 
} // namespace impactx::integrators
 
#endif // IMPACTX_INTEGRATORS_H_