diff options
| author |  Luca Fedeli <luca.fedeli@cea.fr>
| 2021-10-12 00:56:12 +0200 |
|---|---|---|
| committer |  GitHub <noreply@github.com>
| 2021-10-11 15:56:12 -0700 |
| commit | 619fe5e264f9feb192f04d73bf5cfa5fae034abf (patch) | |
| tree | 235eacc5c538b0a974e616c9deefe92fdfbc123d /Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp | |
| parent | 17e157e8eefece6c3e77ab6884fb3f2cd1aa3b8f (diff) | |
| download | WarpX-619fe5e264f9feb192f04d73bf5cfa5fae034abf.tar.gz WarpX-619fe5e264f9feb192f04d73bf5cfa5fae034abf.tar.zst WarpX-619fe5e264f9feb192f04d73bf5cfa5fae034abf.zip | |
Reduced diagnostics: fix calculation of particle energy in single precision (#2397)
* fix kinetic energy calculation in single precision
* really fix the issue...
* add missing _rt
* we don't need m2 now
* remove unused variable
* always perform kinetic energy calculation in double precision
* put w where it was
* fixed bug
* fixed bug
* use threshold for kinetic energy
* fix bug
* increase order of Taylor expansion
* increase order of Taylor expansion
* reduce classical/relativistic threshold for kinetic energy
* reset threshold to 1.005
* improve test
Diffstat (limited to 'Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp')
| -rw-r--r-- | Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp | 16 |
1 files changed, 11 insertions, 5 deletions
diff --git a/Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp b/Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp index b5196a032..7485b5dbf 100644 --- a/Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp +++ b/Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp @@ -101,8 +101,8 @@ void ParticleEnergy::ComputeDiags (int step) const int nSpecies = mypc.nSpecies(); // Some useful constants - amrex::Real c2 = PhysConst::c * PhysConst::c; - amrex::Real c4 = c2 * c2; + constexpr amrex::Real c2 = PhysConst::c * PhysConst::c; + constexpr amrex::Real inv_c2 = 1.0_rt/c2; // Some useful offsets to fill m_data below int offset_total_species, offset_mean_species, offset_mean_all; @@ -117,7 +117,6 @@ void ParticleEnergy::ComputeDiags (int step) // Get mass (used only for particles other than photons, see below) amrex::Real m = myspc.getMass(); - amrex::Real m2 = m * m; using PType = typename WarpXParticleContainer::SuperParticleType; @@ -159,8 +158,15 @@ void ParticleEnergy::ComputeDiags (int step) const amrex::Real ux = p.rdata(PIdx::ux); const amrex::Real uy = p.rdata(PIdx::uy); const amrex::Real uz = p.rdata(PIdx::uz); - const amrex::Real us = ux*ux + uy*uy + uz*uz; - return {w*(std::sqrt(us*m2*c2 + m2*c4) - m*c2), w}; + const amrex::Real u2 = (ux*ux + uy*uy + uz*uz)*inv_c2; + const auto gamma = std::sqrt(1.0_rt + u2); + + const auto kk = (gamma > PhysConst::gamma_relativistic_threshold)? + (gamma-1.0_rt): + (u2*0.5_rt - u2*u2*(1.0_rt/8.0_rt) + u2*u2*u2*(1.0_rt/16.0_rt)- + u2*u2*u2*u2*(5.0_rt/128.0_rt) + (7.0_rt/256_rt)*u2*u2*u2*u2*u2); //Taylor expansion + + return {w*m*c2*kk, w}; }, reduce_ops); |