diff options
| -rw-r--r-- | Docs/source/usage/parameters.rst | 10 | ||||
| -rwxr-xr-x | Examples/Tests/plasma_lens/analysis.py | 102 | ||||
| -rw-r--r-- | Examples/Tests/plasma_lens/inputs_3d | 50 | ||||
| -rw-r--r-- | Regression/Checksum/benchmarks_json/Plasma_lens.json | 23 | ||||
| -rw-r--r-- | Regression/WarpX-tests.ini | 18 | ||||
| -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 |
9 files changed, 300 insertions, 2 deletions
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 99b9e259d..13e6f5319 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -1134,8 +1134,16 @@ Laser initialization using ``my_constants``. For a two-dimensional simulation, similar to the B-field, it is assumed that the first and second dimensions are `x` and `z`, respectively, and the value of the `Ey` component is set to zero. - The current implementation of the parser for B-field on particles + The current implementation of the parser for E-field on particles is applied in cartesian co-ordinates as a function of (x,y,z) even for RZ. + To apply a series of plasma lenses, use the option ``repeated_plasma_lens``. This + option requires the following parameters, + ``repeated_plasma_lens_period``, the period length of the repeat, a single float number, + ``repeated_plasma_lens_starts``, the start of each lens relative to the period, an array of floats, + ``repeated_plasma_lens_lengths``, the length of each lens, an array of floats, + ``repeated_plasma_lens_strengths``, the focusing strength of each lens, an array of floats. + The applied field is uniform longitudinally (along z) with a hard edge, + where residence corrections are used for more accurate field calculation. * ``particles.E_external_particle`` & ``particles.B_external_particle`` (list of `float`) optional (default `0. 0. 0.`) Two separate parameters which add an externally applied uniform E-field or diff --git a/Examples/Tests/plasma_lens/analysis.py b/Examples/Tests/plasma_lens/analysis.py new file mode 100755 index 000000000..b348a3053 --- /dev/null +++ b/Examples/Tests/plasma_lens/analysis.py @@ -0,0 +1,102 @@ +#! /usr/bin/env python + +# Copyright 2021 David Grote +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + + +""" +This script tests the plasma lens. +The input file sets up a series of plasma lens and propagates a particle through them. +The final position is compared to the analytic solution. +The motion is slow enough that relativistic effects are ignored. +""" + +import sys +import os +import yt +import numpy as np +from scipy.constants import e, m_e, c +yt.funcs.mylog.setLevel(0) +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +filename = sys.argv[1] +ds = yt.load( filename ) +ad = ds.all_data() + +# Get final position of the particles +# The sorting by id is to get them in the correct order +# There are two particles, one moves in x, the other in y. +r_id = ad['electrons', 'particle_id'].v + +xx_sim = ad['electrons', 'particle_position_x'].v[np.argsort(r_id)][0] +yy_sim = ad['electrons', 'particle_position_y'].v[np.argsort(r_id)][1] +zz_sim0 = ad['electrons', 'particle_position_z'].v[np.argsort(r_id)][0] +zz_sim1 = ad['electrons', 'particle_position_z'].v[np.argsort(r_id)][1] + +ux_sim = ad['electrons', 'particle_momentum_x'].v[np.argsort(r_id)][0]/m_e +uy_sim = ad['electrons', 'particle_momentum_y'].v[np.argsort(r_id)][1]/m_e + + +def applylens(x0, vx0, vz0, lens_length, lens_strength): + """Given the initial position at the start of the lens, + calculate the end position at the exit. + This assumes harmonic motion.""" + w = np.sqrt(e/m_e*lens_strength) + A = np.sqrt(x0**2 + (vx0/w)**2) + phi = np.arctan2(-vx0, w*x0) + t = lens_length/vz0 + x1 = A*np.cos(w*t + phi) + vx1 = -w*A*np.sin(w*t + phi) + return x1, vx1 + +plasma_lens_period = float(ds.parameters.get('particles.repeated_plasma_lens_period')) +plasma_lens_starts = [float(x) for x in ds.parameters.get('particles.repeated_plasma_lens_starts').split()] +plasma_lens_lengths = [float(x) for x in ds.parameters.get('particles.repeated_plasma_lens_lengths').split()] +plasma_lens_strengths = [eval(x) for x in ds.parameters.get('particles.repeated_plasma_lens_strengths').split()] + +clight = c + +x0 = float(ds.parameters.get('electrons.multiple_particles_pos_x').split()[0]) +y0 = float(ds.parameters.get('electrons.multiple_particles_pos_y').split()[1]) +z0 = float(ds.parameters.get('electrons.multiple_particles_pos_z').split()[0]) +ux0 = float(ds.parameters.get('electrons.multiple_particles_vel_x').split()[0])*c +uy0 = float(ds.parameters.get('electrons.multiple_particles_vel_y').split()[1])*c +uz0 = eval(ds.parameters.get('electrons.multiple_particles_vel_z').split()[0])*c + +tt = 0. +xx = x0 +yy = y0 +zz = z0 +ux = ux0 +uy = uy0 +uz = uz0 + +for i in range(len(plasma_lens_starts)): + z_lens = i*plasma_lens_period + plasma_lens_starts[i] + dt = (z_lens - zz)/uz + tt = tt + dt + xx = xx + dt*ux + yy = yy + dt*uy + xx, ux = applylens(xx, ux, uz, plasma_lens_lengths[i], plasma_lens_strengths[i]) + yy, uy = applylens(yy, uy, uz, plasma_lens_lengths[i], plasma_lens_strengths[i]) + dt = plasma_lens_lengths[i]/uz + tt = tt + dt + zz = z_lens + plasma_lens_lengths[i] + +dt0 = (zz_sim0 - zz)/uz +dt1 = (zz_sim1 - zz)/uz +xx = xx + dt0*ux +yy = yy + dt1*uy + +assert abs(xx - xx_sim) < 0.011, Exception('error in x particle position') +assert abs(yy - yy_sim) < 0.011, Exception('error in y particle position') +assert abs(ux - ux_sim) < 70., Exception('error in x particle velocity') +assert abs(uy - uy_sim) < 70., Exception('error in y particle velocity') + +test_name = os.path.split(os.getcwd())[1] +checksumAPI.evaluate_checksum(test_name, filename) + diff --git a/Examples/Tests/plasma_lens/inputs_3d b/Examples/Tests/plasma_lens/inputs_3d new file mode 100644 index 000000000..fa19685a2 --- /dev/null +++ b/Examples/Tests/plasma_lens/inputs_3d @@ -0,0 +1,50 @@ +# Maximum number of time steps +max_step = 97 + +# number of grid points +amr.n_cell = 16 16 16 + +amr.max_level = 0 + +# Geometry +geometry.coord_sys = 0 # 0: Cartesian +geometry.prob_lo = -1.0 -1.0 0.0 # physical domain +geometry.prob_hi = 1.0 1.0 2.0 + +boundary.field_lo = pec pec pec +boundary.field_hi = pec pec pec +boundary.particle_lo = absorbing absorbing absorbing +boundary.particle_hi = absorbing absorbing absorbing + +warpx.do_pml = 0 +warpx.const_dt = 1.e-6 +warpx.do_electrostatic = labframe + +# Algorithms +algo.particle_shape = 1 + +# particles +particles.species_names = electrons + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "MultipleParticles" +electrons.multiple_particles_pos_x = 0.5 0. +electrons.multiple_particles_pos_y = 0. 0.4 +electrons.multiple_particles_pos_z = 0.05 0.05 +electrons.multiple_particles_vel_x = 0. 0. +electrons.multiple_particles_vel_y = 0. 0. +electrons.multiple_particles_vel_z = 0.02e6/clight 0.02e6/clight +electrons.multiple_particles_weight = 1. 1. + +particles.E_ext_particle_init_style = repeated_plasma_lens +particles.repeated_plasma_lens_period = 0.5 +particles.repeated_plasma_lens_starts = 0.1 0.11 0.12 0.13 +particles.repeated_plasma_lens_lengths = 0.1 0.11 0.12 0.13 +particles.repeated_plasma_lens_strengths = 0.07 0.06 0.06 0.03 + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 97 +diag1.diag_type = Full +diag1.electrons.variables = ux uy uz diff --git a/Regression/Checksum/benchmarks_json/Plasma_lens.json b/Regression/Checksum/benchmarks_json/Plasma_lens.json new file mode 100644 index 000000000..69ad58618 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/Plasma_lens.json @@ -0,0 +1,23 @@ +{ + "electrons": { + "particle_cpu": 1.0, + "particle_id": 2.0, + "particle_momentum_x": 2.2518779586696112e-26, + "particle_momentum_y": 1.8015029214643703e-26, + "particle_momentum_z": 3.6437524496743613e-26, + "particle_position_x": 0.33581536370442583, + "particle_position_y": 0.2686524426674559, + "particle_position_z": 3.979988490885961 + }, + "lev=0": { + "Bx": 0.0, + "By": 0.0, + "Bz": 0.0, + "Ex": 8.557217761445071e-07, + "Ey": 9.086796476605627e-07, + "Ez": 1.2874219492176113e-06, + "jx": 0.0, + "jy": 0.0, + "jz": 0.0 + } +} diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 5ab0c3c0f..f79c7eb8e 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -2188,6 +2188,24 @@ particleTypes = electrons analysisRoutine = Examples/Tests/pass_mpi_communicator/analysis.py tolerance = 1.e-14 +[Plasma_lens] +buildDir = . +inputFile = Examples/Tests/plasma_lens/inputs_3d +runtime_params = +dim = 3 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = electrons +analysisRoutine = Examples/Tests/plasma_lens/analysis.py +tolerance = 1.e-14 + [background_mcc] buildDir = . inputFile = Examples/Physics_applications/capacitive_discharge/inputs_2d 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); |