diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/Particles/Gather/GetExternalFields.H | 51 | ||||
| -rw-r--r-- | Source/Particles/Gather/GetExternalFields.cpp | 15 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.H | 8 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.cpp | 25 |
4 files changed, 98 insertions, 1 deletions
diff --git a/Source/Particles/Gather/GetExternalFields.H b/Source/Particles/Gather/GetExternalFields.H index 254d335b9..468b8d078 100644 --- a/Source/Particles/Gather/GetExternalFields.H +++ b/Source/Particles/Gather/GetExternalFields.H @@ -12,7 +12,7 @@ #include <AMReX_Parser.H> #include <AMReX_REAL.H> -enum ExternalFieldInitType { Constant, Parser }; +enum ExternalFieldInitType { Constant, Parser, RepeatedPlasmaLens }; /** \brief Base class for functors that assign external * field values (E or B) to particles. @@ -29,12 +29,23 @@ struct GetExternalField GetParticlePosition m_get_position; amrex::Real m_time; + amrex::Real m_repeated_plasma_lens_period; + const amrex::Real* AMREX_RESTRICT m_repeated_plasma_lens_starts = nullptr; + const amrex::Real* AMREX_RESTRICT m_repeated_plasma_lens_lengths = nullptr; + const amrex::Real* AMREX_RESTRICT m_repeated_plasma_lens_strengths = nullptr; + int m_n_lenses; + amrex::Real m_dt; + const amrex::ParticleReal* AMREX_RESTRICT m_ux = nullptr; + const amrex::ParticleReal* AMREX_RESTRICT m_uy = nullptr; + const amrex::ParticleReal* AMREX_RESTRICT m_uz = nullptr; + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator () (long i, amrex::ParticleReal& field_x, amrex::ParticleReal& field_y, amrex::ParticleReal& field_z) const noexcept { + using namespace amrex::literals; if (m_type == Constant) { field_x += m_field_value[0]; @@ -49,6 +60,44 @@ struct GetExternalField field_y += m_yfield_partparser(x, y, z, m_time); field_z += m_zfield_partparser(x, y, z, m_time); } + else if (m_type == RepeatedPlasmaLens) + { + amrex::ParticleReal x, y, z; + m_get_position(i, x, y, z); + + amrex::ParticleReal const uxp = m_ux[i]; + amrex::ParticleReal const uyp = m_uy[i]; + amrex::ParticleReal const uzp = m_uz[i]; + constexpr amrex::Real inv_c2 = 1._rt/(PhysConst::c*PhysConst::c); + const amrex::Real inv_gamma = 1._rt/std::sqrt(1._rt + (uxp*uxp + uyp*uyp + uzp*uzp)*inv_c2); + const amrex::ParticleReal vzp = uzp*inv_gamma; + + // This assumes that vzp > 0. + amrex::ParticleReal const zl = z; + amrex::ParticleReal const zr = z + vzp*m_dt; + + // This assumes that zl > 0. + int i_lens = static_cast<int>(std::floor(zl/m_repeated_plasma_lens_period)); + i_lens = i_lens % m_n_lenses; + amrex::Real const lens_start = m_repeated_plasma_lens_starts[i_lens] + i_lens*m_repeated_plasma_lens_period; + amrex::Real const lens_end = lens_start + m_repeated_plasma_lens_lengths[i_lens]; + + // Calculate the residence correction + // frac will be 1 if the step is completely inside the lens, between 0 and 1 + // when entering or leaving the lens, and otherwise 0. + amrex::Real fl = 0.; + if (zl >= lens_start && zl < lens_end) fl = 1.; + amrex::Real fr = 0.; + if (zr >= lens_start && zr < lens_end) fr = 1.; + amrex::Real frac = fl; + amrex::Real dzi = 1./(vzp*m_dt); + if (fl > fr) frac = (lens_end - zl)*dzi; + if (fr > fl) frac = (zr - lens_start)*dzi; + + field_x += x*frac*m_repeated_plasma_lens_strengths[i_lens]; + field_y += y*frac*m_repeated_plasma_lens_strengths[i_lens]; + + } else { amrex::Abort("ExternalFieldInitType not known!!! \n"); diff --git a/Source/Particles/Gather/GetExternalFields.cpp b/Source/Particles/Gather/GetExternalFields.cpp index 0808bd6c1..03c4dc2b8 100644 --- a/Source/Particles/Gather/GetExternalFields.cpp +++ b/Source/Particles/Gather/GetExternalFields.cpp @@ -29,6 +29,21 @@ GetExternalEField::GetExternalEField (const WarpXParIter& a_pti, int a_offset) n m_yfield_partparser = mypc.m_Ey_particle_parser->compile<4>(); m_zfield_partparser = mypc.m_Ez_particle_parser->compile<4>(); } + else if (mypc.m_E_ext_particle_s=="repeated_plasma_lens") + { + m_type = RepeatedPlasmaLens; + m_dt = warpx.getdt(a_pti.GetLevel()); + m_get_position = GetParticlePosition(a_pti, a_offset); + auto& attribs = a_pti.GetAttribs(); + m_ux = attribs[PIdx::ux].dataPtr() + a_offset; + m_uy = attribs[PIdx::uy].dataPtr() + a_offset; + m_uz = attribs[PIdx::uz].dataPtr() + a_offset; + m_repeated_plasma_lens_period = mypc.m_repeated_plasma_lens_period; + m_n_lenses = static_cast<int>(mypc.h_repeated_plasma_lens_starts.size()); + m_repeated_plasma_lens_starts = mypc.d_repeated_plasma_lens_starts.data(); + m_repeated_plasma_lens_lengths = mypc.d_repeated_plasma_lens_lengths.data(); + m_repeated_plasma_lens_strengths = mypc.d_repeated_plasma_lens_strengths.data(); + } } GetExternalBField::GetExternalBField (const WarpXParIter& a_pti, int a_offset) noexcept diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 5ee416add..b59aa4371 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -296,6 +296,14 @@ public: std::unique_ptr<amrex::Parser> m_Ey_particle_parser; std::unique_ptr<amrex::Parser> m_Ez_particle_parser; + amrex::Real m_repeated_plasma_lens_period; + amrex::Vector<amrex::Real> h_repeated_plasma_lens_starts; + amrex::Vector<amrex::Real> h_repeated_plasma_lens_lengths; + amrex::Vector<amrex::Real> h_repeated_plasma_lens_strengths; + amrex::Gpu::DeviceVector<amrex::Real> d_repeated_plasma_lens_starts; + amrex::Gpu::DeviceVector<amrex::Real> d_repeated_plasma_lens_lengths; + amrex::Gpu::DeviceVector<amrex::Real> d_repeated_plasma_lens_strengths; + #ifdef WARPX_QED /** * \brief Performs QED events (Breit-Wheeler process and photon emission) diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 5e59d6eaa..b0a1d211f 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -223,6 +223,31 @@ MultiParticleContainer::ReadParameters () } + // if the input string for E_ext_particle_s is + // "repeated_plasma_lens" then the plasma lens properties + // must be provided in the input file. + if (m_E_ext_particle_s == "repeated_plasma_lens") { + queryWithParser(pp_particles, "repeated_plasma_lens_period", m_repeated_plasma_lens_period); + getArrWithParser(pp_particles, "repeated_plasma_lens_starts", h_repeated_plasma_lens_starts); + getArrWithParser(pp_particles, "repeated_plasma_lens_lengths", h_repeated_plasma_lens_lengths); + getArrWithParser(pp_particles, "repeated_plasma_lens_strengths", h_repeated_plasma_lens_strengths); + + int n_lenses = static_cast<int>(h_repeated_plasma_lens_starts.size()); + d_repeated_plasma_lens_starts.resize(n_lenses); + d_repeated_plasma_lens_lengths.resize(n_lenses); + d_repeated_plasma_lens_strengths.resize(n_lenses); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + h_repeated_plasma_lens_starts.begin(), h_repeated_plasma_lens_starts.end(), + d_repeated_plasma_lens_starts.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + h_repeated_plasma_lens_lengths.begin(), h_repeated_plasma_lens_lengths.end(), + d_repeated_plasma_lens_lengths.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + h_repeated_plasma_lens_strengths.begin(), h_repeated_plasma_lens_strengths.end(), + d_repeated_plasma_lens_strengths.begin()); + amrex::Gpu::synchronize(); + } + // particle species pp_particles.queryarr("species_names", species_names); |