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
|
/* Copyright 2019 David Grote, Maxence Thevenet, Remi Lehe
* Weiqun Zhang
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#ifndef WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_VAY_H_
#define WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_VAY_H_
#include "Utils/WarpXConst.H"
#include <AMReX_FArrayBox.H>
#include <AMReX_REAL.H>
/** \brief Push the particle's positions over one timestep,
* given the value of its momenta `ux`, `uy`, `uz` */
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void UpdateMomentumVay(
amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz,
const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez,
const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz,
const amrex::ParticleReal q, const amrex::ParticleReal m, const amrex::Real dt )
{
using namespace amrex::literals;
// Constants
const amrex::ParticleReal econst = q*dt/m;
const amrex::ParticleReal bconst = 0.5_prt*q*dt/m;
constexpr amrex::ParticleReal invclight = 1._prt/PhysConst::c;
constexpr amrex::ParticleReal invclightsq = 1._prt/(PhysConst::c*PhysConst::c);
// Compute initial gamma
const amrex::ParticleReal inv_gamma = 1._prt/std::sqrt(1._prt + (ux*ux + uy*uy + uz*uz)*invclightsq);
// Get tau
const amrex::ParticleReal taux = bconst*Bx;
const amrex::ParticleReal tauy = bconst*By;
const amrex::ParticleReal tauz = bconst*Bz;
const amrex::ParticleReal tausq = taux*taux+tauy*tauy+tauz*tauz;
// Get U', gamma'^2
const amrex::ParticleReal uxpr = ux + econst*Ex + (uy*tauz-uz*tauy)*inv_gamma;
const amrex::ParticleReal uypr = uy + econst*Ey + (uz*taux-ux*tauz)*inv_gamma;
const amrex::ParticleReal uzpr = uz + econst*Ez + (ux*tauy-uy*taux)*inv_gamma;
const amrex::ParticleReal gprsq = (1._prt + (uxpr*uxpr + uypr*uypr + uzpr*uzpr)*invclightsq);
// Get u*
const amrex::ParticleReal ust = (uxpr*taux + uypr*tauy + uzpr*tauz)*invclight;
// Get new gamma
const amrex::ParticleReal sigma = gprsq-tausq;
const amrex::ParticleReal gisq = 2._prt/(sigma + std::sqrt(sigma*sigma + 4._prt*(tausq + ust*ust)) );
// Get t, s
const amrex::ParticleReal bg = bconst*std::sqrt(gisq);
const amrex::ParticleReal tx = bg*Bx;
const amrex::ParticleReal ty = bg*By;
const amrex::ParticleReal tz = bg*Bz;
const amrex::ParticleReal s = 1._prt/(1._prt+tausq*gisq);
// Get t.u'
const amrex::ParticleReal tu = tx*uxpr + ty*uypr + tz*uzpr;
// Get new U
ux = s*(uxpr+tx*tu+uypr*tz-uzpr*ty);
uy = s*(uypr+ty*tu+uzpr*tx-uxpr*tz);
uz = s*(uzpr+tz*tu+uxpr*ty-uypr*tx);
}
#endif // WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_VAY_H_
|
