BorisPusher.h 3.48 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
#ifndef CLASSIC_PartPusher_H
#define CLASSIC_PartPusher_H

#include "Algorithms/Vektor.h"
#include "Algorithms/PartData.h"
#include "Physics/Physics.h"

/*


 */

class BorisPusher {
    
public:
    BorisPusher(const PartData &ref);
    BorisPusher();
    void initialise(const PartData *ref);
    
    void kick(const Vector_t &R, Vector_t &P,
              const Vector_t &Ef, const Vector_t &Bf,
              const double &dt) const;
    
    void kick(const Vector_t &R, Vector_t &P,
              const Vector_t &Ef, const Vector_t &Bf,
              const double &dt, const double &mass,
              const double &charge) const;
    
    void push(Vector_t &R, const Vector_t &P, const double &dt) const;

private:
    const PartData *itsReference;
};

inline BorisPusher::BorisPusher(const PartData &ref):
    itsReference(&ref)
{ }

inline BorisPusher::BorisPusher():
    itsReference(NULL)
{ }

inline void BorisPusher::initialise(const PartData *ref)
{
    itsReference = ref;
}

inline void BorisPusher::kick(const Vector_t &R, Vector_t &P,
                              const Vector_t &Ef, const Vector_t &Bf,
                              const double &dt) const
{
    kick(R, P, Ef, Bf, dt, itsReference->getM(), itsReference->getQ());
}


inline void BorisPusher::kick(const Vector_t &R, Vector_t &P,
                              const Vector_t &Ef, const Vector_t &Bf,
                              const double &dt, const double &mass,
                              const double &charge) const
{
    // Implementation follows chapter 4-4, p. 61 - 63 from
    // Birdsall, C. K. and Langdon, A. B. (1985). Plasma physics
    // via computer simulation.
    //
    // Up to finite precision effects, the new implementation is equivalent to the
    // old one, but uses less floating point operations.
    //
    // Relativistic variant implemented below is described in
    // chapter 15-4, p. 356 - 357.
    // However, since other units are used here, small
    // modifications are required. The relativistic variant can be derived
    // from the nonrelativistic one by replacing
    //     mass
    // by
    //     gamma * rest mass
    // and transforming the units.
    //
    // Parameters:
    //     R = x / (c * dt): Scaled position x, not used in here
    //     P = v / c * gamma: Scaled velocity v
    //     Ef: Electric field
    //     Bf: Magnetic field
    //     dt: Timestep
    //     mass = rest energy = rest mass * c * c
    //     charge

    using Physics::c;

    // Half step E
    P += 0.5 * dt * charge * c / mass * Ef;

    // Full step B

    /*
        LF
    double const gamma = sqrt(1.0 + dot(P, P));
    Vector_t const t = dt * charge * c * c / (gamma * mass) * Bf;
    P +=  cross(P, t);
    */

    double const gamma = sqrt(1.0 + dot(P, P));
    Vector_t const t = 0.5 * dt * charge * c * c / (gamma * mass) * Bf;
    Vector_t const w = P + cross(P, t);
    Vector_t const s = 2.0 / (1.0 + dot(t, t)) * t;
    P += cross(w, s);

    /* a poor Leap-Frog
        P += 1.0 * dt * charge * c * c / (gamma * mass) * cross(P,Bf);
    */

    // Half step E
    P += 0.5 * dt * charge * c / mass * Ef;
}


inline void BorisPusher::push(Vector_t &R, const Vector_t &P, const double &/* dt */) const {
    /** \f[ \vec{x}_{n+1/2} = \vec{x}_{n} + \frac{1}{2}\vec{v}_{n-1/2}\quad (= \vec{x}_{n} + \frac{\Delta t}{2} \frac{\vec{\beta}_{n-1/2}\gamma_{n-1/2}}{\gamma_{n-1/2}}) \f]
     *
     * \code
     * R[i] += 0.5 * P[i] * recpgamma;
     * \endcode
     */
    R += 0.5 * P / sqrt(1.0 + dot(P, P));
}

#endif