Added B field to plasma lens (#2163)

* Added B field to plasma lens

* Fix B field and updated CI test to include the B field

* Updated benchmark

Examples/Tests/plasma_lens/analysis.py

@@ -53,10 +53,13 @@ def applylens(x0, vx0, vz0, lens_length, lens_strength):
     vx1 = -w*A*np.sin(w*t + phi)
     return x1, vx1
 
+vel_z = eval(ds.parameters.get('my_constants.vel_z'))
+
 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()]
+plasma_lens_strengths_E = [eval(x) for x in ds.parameters.get('particles.repeated_plasma_lens_strengths_E').split()]
+plasma_lens_strengths_B = [eval(x) for x in ds.parameters.get('particles.repeated_plasma_lens_strengths_B').split()]
 
 clight = c
 
@@ -81,8 +84,9 @@ def applylens(x0, vx0, vz0, lens_length, lens_strength):
     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])
+    lens_strength = plasma_lens_strengths_E[i] + plasma_lens_strengths_B[i]*vel_z
+    xx, ux = applylens(xx, ux, uz, plasma_lens_lengths[i], lens_strength)
+    yy, uy = applylens(yy, uy, uz, plasma_lens_lengths[i], lens_strength)
     dt = plasma_lens_lengths[i]/uz
     tt = tt + dt
     zz = z_lens + plasma_lens_lengths[i]
@@ -92,10 +96,10 @@ def applylens(x0, vx0, vz0, lens_length, lens_strength):
 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')
+assert abs(np.abs((xx - xx_sim)/xx)) < 0.003, Exception('error in x particle position')
+assert abs(np.abs((yy - yy_sim)/yy)) < 0.003, Exception('error in y particle position')
+assert abs(np.abs((ux - ux_sim)/ux)) < 5.e-5, Exception('error in x particle velocity')
+assert abs(np.abs((uy - uy_sim)/uy)) < 5.e-5, Exception('error in y particle velocity')
 
 test_name = os.path.split(os.getcwd())[1]
 checksumAPI.evaluate_checksum(test_name, filename)

Examples/Tests/plasma_lens/inputs_3d

@@ -17,31 +17,35 @@ boundary.particle_lo = absorbing absorbing absorbing
 boundary.particle_hi = absorbing absorbing absorbing
 
 warpx.do_pml = 0
-warpx.const_dt = 1.e-6
+warpx.const_dt = 1.e-7
 warpx.do_electrostatic = labframe
 
 # Algorithms
 algo.particle_shape = 1
 
+my_constants.vel_z = 0.2e6
+
 # 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_x = 0.05 0.
+electrons.multiple_particles_pos_y = 0. 0.04
 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_vel_z = vel_z/clight vel_z/clight
 electrons.multiple_particles_weight = 1. 1.
 
 particles.E_ext_particle_init_style = repeated_plasma_lens
+particles.B_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
+particles.repeated_plasma_lens_strengths_E = 0.06 0.08 0.06 0.02
+particles.repeated_plasma_lens_strengths_B = 0.08/vel_z 0.04/vel_z 0.06/vel_z 0.04/vel_z
 
 # Diagnostics
 diagnostics.diags_names = diag1

Regression/Checksum/benchmarks_json/Plasma_lens.json

@@ -2,20 +2,20 @@
   "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
+    "particle_momentum_x": 1.8497357687843446e-27,
+    "particle_momentum_y": 1.4797852897079812e-27,
+    "particle_momentum_z": 3.6436751241117655e-25,
+    "particle_position_x": 0.03778437270526418,
+    "particle_position_y": 0.030227636257933975,
+    "particle_position_z": 3.9799632585096267
   },
   "lev=0": {
     "Bx": 0.0,
     "By": 0.0,
     "Bz": 0.0,
-    "Ex": 8.557217761445071e-07,
-    "Ey": 9.086796476605627e-07,
-    "Ez": 1.2874219492176113e-06,
+    "Ex": 9.938634536083956e-07,
+    "Ey": 9.994968813016262e-07,
+    "Ez": 1.3536427852197406e-06,
     "jx": 0.0,
     "jy": 0.0,
     "jz": 0.0

Source/Particles/Gather/GetExternalFields.H

@@ -34,6 +34,7 @@ struct GetExternalField
     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;
+    int m_lens_is_electric;
     amrex::Real m_dt;
     const amrex::ParticleReal* AMREX_RESTRICT m_ux = nullptr;
     const amrex::ParticleReal* AMREX_RESTRICT m_uy = nullptr;
@@ -94,8 +95,13 @@ struct GetExternalField
             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];
+            if (m_lens_is_electric) {
+                field_x += x*frac*m_repeated_plasma_lens_strengths[i_lens];
+                field_y += y*frac*m_repeated_plasma_lens_strengths[i_lens];
+            } else {
+                field_x += y*frac*m_repeated_plasma_lens_strengths[i_lens];
+                field_y -= x*frac*m_repeated_plasma_lens_strengths[i_lens];
+            }
 
         }
         else

Source/Particles/Gather/GetExternalFields.cpp

@@ -32,6 +32,7 @@ GetExternalEField::GetExternalEField (const WarpXParIter& a_pti, int a_offset) n
     else if (mypc.m_E_ext_particle_s=="repeated_plasma_lens")
     {
         m_type = RepeatedPlasmaLens;
+        m_lens_is_electric = 1;
         m_dt = warpx.getdt(a_pti.GetLevel());
         m_get_position = GetParticlePosition(a_pti, a_offset);
         auto& attribs = a_pti.GetAttribs();
@@ -42,7 +43,7 @@ GetExternalEField::GetExternalEField (const WarpXParIter& a_pti, int a_offset) n
         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();
+        m_repeated_plasma_lens_strengths = mypc.d_repeated_plasma_lens_strengths_E.data();
     }
 }
 
@@ -66,4 +67,20 @@ GetExternalBField::GetExternalBField (const WarpXParIter& a_pti, int a_offset) n
         m_yfield_partparser = mypc.m_By_particle_parser->compile<4>();
         m_zfield_partparser = mypc.m_Bz_particle_parser->compile<4>();
     }
+    else if (mypc.m_B_ext_particle_s=="repeated_plasma_lens")
+    {
+        m_type = RepeatedPlasmaLens;
+        m_lens_is_electric = 0;
+        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_B.data();
+    }
 }

Source/Particles/MultiParticleContainer.H

@@ -299,10 +299,12 @@ public:
     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::Vector<amrex::Real> h_repeated_plasma_lens_strengths_E;
+    amrex::Vector<amrex::Real> h_repeated_plasma_lens_strengths_B;
     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;
+    amrex::Gpu::DeviceVector<amrex::Real> d_repeated_plasma_lens_strengths_E;
+    amrex::Gpu::DeviceVector<amrex::Real> d_repeated_plasma_lens_strengths_B;
 
 #ifdef WARPX_QED
     /**

Source/Particles/MultiParticleContainer.cpp

@@ -223,28 +223,39 @@ MultiParticleContainer::ReadParameters ()
 
         }
 
-        // if the input string for E_ext_particle_s is
+        // if the input string for E_ext_particle_s or B_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") {
+        if (m_E_ext_particle_s == "repeated_plasma_lens" ||
+            m_B_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());
+
+            if (m_E_ext_particle_s == "repeated_plasma_lens") {
+                getArrWithParser(pp_particles, "repeated_plasma_lens_strengths_E", h_repeated_plasma_lens_strengths_E);
+                d_repeated_plasma_lens_strengths_E.resize(n_lenses);
+                amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+                           h_repeated_plasma_lens_strengths_E.begin(), h_repeated_plasma_lens_strengths_E.end(),
+                           d_repeated_plasma_lens_strengths_E.begin());
+            }
+            if (m_B_ext_particle_s == "repeated_plasma_lens") {
+                getArrWithParser(pp_particles, "repeated_plasma_lens_strengths_B", h_repeated_plasma_lens_strengths_B);
+                d_repeated_plasma_lens_strengths_B.resize(n_lenses);
+                amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+                           h_repeated_plasma_lens_strengths_B.begin(), h_repeated_plasma_lens_strengths_B.end(),
+                           d_repeated_plasma_lens_strengths_B.begin());
+            }
             amrex::Gpu::synchronize();
         }